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1  SUMMARY 


This  report  presents  the  results  of  an  experimental  investigation  into  the  combined 
effects  of  inelasticity  and  strain  rate  sensitivity  on  penetration  into  geologic  or  geo¬ 
logical  derived  targets.  Both  material  models  and  specific  computational  methods 
have  been  developed.  A  new  class  of  compressible  rigid  viscoplastic  models  were  pro¬ 
posed  to  capture  the  solid-fluid  transition  in  behavior  at  high  strain  rates,  account 
for  damage/ plasticity  couplings  and  viscous  effects  which  are  observed  in  geological 
and  cementitious  materials, 

A  hybrid  time-discretization  was  used  to  model  the  non-station  ary  flow  of  the  the 
target  material  and  the  projectile- target  interaction,  i.c,  an  explicit  Euler  method 
for  the  projectile  equation  and  a  forward  (implicit)  method  for  the  target  bound¬ 
ary  value  problem.  At  each  time  step,  a  mixed  finite-element  and  finite- volume 
strategy  was  used  to  solve  the  ” target”  boundary  value  problem.  Specifically,  the 
nonlinear  variational  inequality  for  the  velocity  field  was  discretized  using  the  fi¬ 
nite  element  method  while  a  finite  volume  method  was  used  for  the  hyperbolic 
mass  conservation  and  damage  evolution  equations.  To  solve  the  velocity  problem, 
a  decomposition-coordination  formulation  coupled  with  the  augmented  lagrangian 
method  was  adopted. 

Numerical  simulations  of  penetration  into  concrete  were  performed.  By  conduct¬ 
ing  a  time  step  sensitivity  study  it  was  shown  that  the  numerical  model  is  robust  and 
computationally  inexpensive.  For  the  constants  involved  in  the  model  (shear  and 
volumetric  viscosities,  cut-off  yield  limit  and  exponential  weakening  parameter  for 
friction)  that  cannot  be  determined  from  data,  a  parametric  study  was  performed. 
It  is  shown  that  using  the  material  model  and  numerical  algorithms  developed  the 
evolution  of  the  density  changes  around  the  penetration  tunnel,  the  shape  and  lo¬ 
cation  of  the  rigid/plastic  boundary,  of  the  compaction  zones  and  of  the  extent  of 
damage  due  to  air  void  collapse  arc  described  accurately. 
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2  Task  1:  Development  of  a  new  class  of  fluid-like 
constitutive  models  for  porous  materials 

2.1  Introduction 

A  general  framework  for  modeling  the  non-linear  behavior  of  cementitious  or  geolog¬ 
ical  materials  for  applications  involving  dynamic  loading  is  offered  by  viscoplasticity 
theory.  In  the  clastic/ viscoplastic  models  with  non-station  ary  yield  surface  the  cen¬ 
tral  idea  is  that  the  loading  function  depends  on  the  stress  state,  the  plastic  history 
of  the  material  through  the  plastic  strains,  and  a  scalar  parameter  which  embodies 
time  dependent  alteration  of  the  material  properties.  A  very  large  number  of  models 
developed  for  inviscld  behavior  have  been  extended  into  the  rate-dependent  range 
using  the  overstress  concept.  However,  most  of  these  models  have  been  developed 
in  support  of  civil  engineering  applications  and  thus  are  in  general  representative  of 
relatively  small  strains  and  low  strain  rate  behavior. 

In  some  processes,  such  as  impact  or  high  velocity  penetration  into  cementitious 
or  geological  materials,  large  deformations,  high  pressures  and  high  strain  rates  oc¬ 
cur.  That  is  why  in  the  mechanical  and  numerical  modeling,  mi  eulcrian  description 
and  a  fluid-type  constitutive  law  seems  to  be  more  suitable.  Fluid-type  constitu¬ 
tive  equations  have  been  used  to  describe  the  high-strain  rate  behavior  of  metallic 
materials  (sec  for  instance  Batra,  1987;  Batra  and  Li,  1988).  Implicit  in  these  mod¬ 
els  is  the  hypothesis  of  incompressibility,  hence  those  models  cannot  describe  the 
observed  irreversible  volumetric  changes  observed  in  geological  and  cementitious  ma¬ 
terials.  To  account  for  the  effect  of  compaction  on  the  response  of  cementitious  or 
geological  materials,  an  explicit  dependence  of  both  the  yield  limit  and  the  viscosity 
coefficients  cm  the  current  density  has  to  be  considered  (  e.g.  Cazacu  ct  ah,  2004). 

It  is  well  recognized  that  depending  on  their  dynamical  state,  particulate  ma¬ 
terials  may  display  both  solid  and  fluid  like  behavior.  Since  at  low  deformation 
rates  the  behavior  is  mainly  dictated  by  frictional  interpart ide  contact  and  sliding 
of  the  particles  over  each  other,  a  large  body  of  literature  exists  in  which  particulate 
systems  arc  modeled  using  a  solid  mechanics  approach.  These  models  arc  gener¬ 
ally  developed  within  the  framework  or  plasticity  theory.  Again,  in  order  to  model 
short-term  time  effects  on  the  quasi-static  deformation  of  such  systems,  extensions 
of  the  inviscid  models  using  Pcrzyna’s  approach  (Perzyna,1966)  have  been  proposed. 
However,  most  of  these  models  are  applicable  only  for  relatively  small  strains  and 
low  strain  rates.  On  the  other  hand,  rapid  unfluidized  flows  such  as  chute  flow's 
and  avalanches  in  which  particles  interact  by  short  collisions  and  cascade  over  each 
other,  are  generally  modeled  based  on  analogy  to  gas  dynamics  {c.g,  Jenkins  and 
Askari,  1991;  Campbell  and  Brennen,  1985)* 

Non-linear  viscous  compressible  equations  have  been  also  used  for  prediction  of 
chute  flows  as  well  as  for  describing  slowly- moving  natural  slopes  which  generally 
have  a  displacement  history  of  few  millimeters  per  year  followed  by  periods  of  ac¬ 
celeration  degenerating  into  catastrophic  events  (e*g,  Viilliet,  2000),  A  combination 
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of  the  two  approaches  was  also  proposed  (e.g.  Savage,  1979;  Goodman  and  Cowin 
(1971).  For  example,  the  premise  of  the  model  proposed  by  Goodman  and  Cowin 
(1971)  is  that  at  rest  the  granular  material  obeys  a  Coulomb  yield  condition.  Once 
the  flow  develops  the  state  of  stress  is  a  function  of  the  rate  of  deformation  and  the 
change  in  volume  distribution,  with  the  dissipative  part  of  the  stress  (i.e,  the  part  of 
stress  that  depends  on  the  deformation  rate)  being  that  of  a  viscous  fluid.  Star  ting 
from  the  non-linear  viscous  model  of  Savage(1979),  recently  Elaskar  and  Godoy 
(1998)  have  proposed  a  non-linear  compressible  viscous  model  of  granular  flow  with 
second-order  effects. 

Experimental  evidence  indicates  that  for  a  wide  range  of  porous  and  particulate 
materials  flow  develops  only  if  a  certain  threshold  that  depends  on  the  stress  level 
and  the  density  or  the  solids  volume  fraction  is  exceeded.  At  difference  with  a 
classical  fluid  model T  in  the  absence  of  motion  such  materials  can  sustain  a  shear 
stress.  Generally,  in  order  to  describe  the  behavior  of  such  materials,  it  is  necessary 
to  formulate  a  flow  activation  condition  or  yield  condition  that  describes  the  state 
of  stress  in  the  material  at  the  onset  of  the  deformation  process  and  to  describe  the 
state  of  stress  in  the  material,  once  the  flow/deformation  develops.  Moreover,  the 
elastic  effects  can  be  neglected, hypothesis  which  is  usually  made  when  describing 
the  large  deformation  behavior  of  particulate  materials  or  cementitious  materials. 

The  main  objective  of  this  work  is  to  provide  procedures  for  constructing  com¬ 
pressible  viscoplastic  fluid- type  constitutive  models  compatible  with  a  given  flow 
initiation  criterion.  Two  methods  will  be  considered:  Perzyna^s  viscoplastic  regu¬ 
larization,  and  a  Maxwell  type  stress  superposition  method.  Examples  of  consti¬ 
tutive  equations  obtained  from  classical  yield  functions  using  both  procedures  are 
presented.  Illustrations  of  the  response  of  these  models  for  uniaxial  compression 
conditions  as  well  as  comparisons  between  these  various  models  arc  provided.  Fur¬ 
thermore,  we  show  that  irrespective  of  the  method  used,  for  certain  yield  condit  ions 
the  resulting  mechanical  model  is  physically  unacceptable.  That  means  that  not 
from  all  yield  conditions  can  a  fluid-type  model  be  deduced.  Extensions  of  these 
models  that  include  coupling  between  damage,  plastic  and  viscous  effects  are  pro¬ 
posed.  Some  results  of  numerical  modeling  of  penetration  into  concrete  using  a 
compressible  rigid- viscoplastic  model  are  presented.  Finally,  extensions  of  fluid- 
type  models  that  include  coupling  between  damage,  plastic  and  viscous  effects  arc 
proposed. 

2.2  Rigid  viscoplastic  fluids 

Let  us  begin  by  describing  what  wc  mean  by  a  rigid  viscoplastic  fluid.  In  contrast 
with  a  classical  fluid,  which  cannot  sustain  a  shear  stress,  we  suppose  that  at  rest 
the  Cauchy  stress  tensor  a  must  belong  to  an  admissible  convex  set  A'  of  R|*3,  the 
space  of  second  order  symmetric  tensors  over  the  set  of  real  numbers  R.  Conversely, 
if  the  stress  is  in  A'  then  the  rate  of  deformation  tensor  D  =  D(v)  =  (Vr  +  V/  v)/2 
(tj  denotes  the  cuierian  velocity  field)  vanishes.  If  the  stress  tensor  is  not  in  A"  then 
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there  isi  How,  the  rate  of  deformation  tensor  D  being  subject  to  certain  kinematic 
constraints,  i.e.  D  €  C,  and  the  stress  tensor  is  a  function  of  D 
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D(v)  e  e, 


tr  =  /(/?(«)) 
(t  e  K 


if|£»(w)|?io, 
if  |X>(tt)|  =0, 


(1) 


It  should  be  noted  that  at  difference  with  a  classical  fluid  constitutive  equation, 
for  a  rigid  viscoplastic  material  the  function  /  is  not  defined  and  cannot  be  prolonged 
by  continuity  in  D  =  0.  To  ensure  continuous  transition  between  flow  and  no-flow 
states,  it  is  necessary  to  impose  that  only  the  limit  points  of  f(D)  when  D  — ►  0  arc 
in  K\  i.e. 

if  D  — *  0  with  D  €  C  and  f[D)  —+  a  then  a  €  K.  (2) 

The  condition  (2)  implies  that  we  neglect  the  elastic  properties  of  the  material.  Such 
a  hypothesis  is  valid  lit  the  high  pressure  regime  in  which  elastic  strains  art?  much 
smaller  than  inelastic  strains.  Also,  what  we  usually  call  the  ” static  pressure"  (  i.e. 
the  pressure  that  the  material  may  sustain  for  zero  volume  deformation  rate)  is  not 
directly  present  in  the  model.  In  order  to  capture  the  effect  of  the  porosity  on  the 
behavior  of  geological  materials  (  c.g.  pore  closure,  pore  collapse),  the  convex  sei  K 
and  /  are  considered  to  be  functions  of  h ,  an  internal  state  variable  related  to  the 
actual  irreversible  volume  change.  Since  elastic  effects  arc  neglected,  dh  —  div{v)dt 
and  based  on  the  continuity  equation 

‘di  =  lH  +  v ' V/>  =  ~pdiv^  (;i) 

where  p  is  the  actual  density,  it  follows  that  d ln(p)  ~dh.  Thus,  the  actual  density 
can  be  considered  as  internal  variable  (see  the  last  section  for  other  choices  of  the 
internal  state  variable  h ),  In  the  following,  the  convex  A"  —  K(p)  will  be  defined 
through  a  continuous  scalar  function  F  =  F(o-,p)  that  describes  the  flow /no  flow 
condition,  i.e, 

K(p)  ={cr€  ;  F(<r.p)  <  0}. 

From  classic  representation  theorems  for  isotropic  tensor  functions,  it  follows 
that  there  exist  three  scalar  functions  a,/3, 7  which  depend  on  D  through  its  three 
invariants  such  that 


f(D,p)  =  a(D,p)I  +  p(D,p)D'  +  7  (D,p)D2. 

where  /  is  the  second-order  identity  tensor,  while  Af  denotes  the  deviator  of  any 
symmetric  second-order  tensor  A ,  i.e.  A*  =  A  —  fr(A)/3J,  and  * trk  stands  for  the 
trace  operator. 

In  the  following,  we  will  neglect  second  order  effects,  i.e.  impose  7  =  0  and 
consider  that  the  constitutive  functions  do  not  depend  on  the  third  invariant  of  the 
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rate  of  deformation  tensor,  Henceforth,  our  discussion  will  be  limited  to  a  subclass 
of  viscoplastic  fluids  characterized  by: 

<r  =  a(D(v),p)I  +  0(D{v),p)D'(v)  if  |D(v)|  #  0, 
F(tr.p)<  0  if  |Z>(«)|  —  0, 

where  C  C  R|x,t  represents  a  kinematic  constraint. 

The  above  constitutive  equation  completely  describes  a  rigid  viscoplastic  fluid 
through  only  four  elements  ;  a  (scalar)  yield  function  F,  a  set  C  of  admissible  strain 
rates  ami  two  scalar  functions  a  and  0.  As  expected,  these  four  elements  are  not 
independent.  In  order  to  define  a  consistent  rigid  viscoplastic  model  compatibility 
conditions  ought  to  be  imposed. 

The  first  two  compatibility  conditions  concern  the  How  regime.  Thus,  if  there 
is  flow  the  associated  stress  a  has  to  be  outside  the  convex  K  defined  through  the 
yield  function  F.  In  this  case  (i.e.  in  the  viscoplastic  regime  where  D  ^  (1)  we  have 
to  impose  that 


D(v)  G  C 


F(a(Dip)I +  0(D,p)D\p)>  0,  for  all  DeC\  0.  (5) 

This  assumption  is  natural  since  the  elastic  effects  arc  neglected  (i.e.  a  rigid  subdo¬ 
main  is  included  in  the  constitutive  domain). 

The  second  compatibility  condition  seems  also  to  be  natural.  Since  wc  model  a 
fluid-type  behavior,  we  suppose  that  a  is  continuous  with  respect  to  D ,  he* 

D  — ►  a(T>,  p),  D  — ►  0(D.  p)  are  continuous  on  C  \  0*  (6) 

That  means  that  small  variations  in  the  rate  of  deformation  D  ^  0  will  give  small 
variations  in  the  stress  field.  Let  us  stress  again  that  for  a  rigid- viscoplastic  fluid 
this  property  is  not  valid  in  the  neighborhood  of  D  =  0. 

The  third  condition  deals  with  the  "continuity”  of  the  transition  between  flow 
and  no  How'.  More  precisely,  the  stress  associated  to  a  fluid  in  motion  which  tends 
to  stop  flowing  has  to  approach  the  yield  surface  (Le,  the  boundary  of  K).  For  the 
particular  form  of  the  constitutive  equation  (65),  the  condition  (2)  becomes 

F(a(D,  p)l  4-  f3(D,p)D\ p)  —  (1,  for  D  0  with  D  €  C  \  0.  (7) 

The  fourth  and  last  compatibility  condition  yields  from  thermodynamic  argu¬ 
ments*  It  states  that  the  dissipated  power  during  viscoplastic  flow  is  positive 
tr  :  D  >  0,  i.e. 


a[D,p)tr{D)  +  0{D.  p)\D'\2  >  0  for  all  D.  (8) 

Our  main  goal  is  to  provide  procedures  for  determining  the  expressions  of  the 
constitutive  functions  F,  ft,  and  0  such  that  the  model  (65)  is  consistent,  (i.e.  the 
compatibility  conditions  (5-8)  arc  satisfied).  Specifically,  we  will  use  two  methods. 
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The  first,  one  is  Perzyna’s  visco-  plastic  regularization  approach  while  the  second 
one  is  based  on  a  Maxwell  type  superposition  of  a  rigid  plastic  material  on  a  com¬ 
pressible  viscous  fluid.  We  illustrate  the  proposed  approach  by  providing  examples 
of  constitutive  equations  obtained  starting  from  classic  yield  functions  by  making 
use  of  both  procedures.  In  general >  even  if  wc  start  from  the  same  yield  condition, 
the  models  obtained  using  these  two  methods  are  different. 

Furthermore,  wc  show  that  starting  from  certain  classical  yield  functions,  we 
can  construct  consistent  fluid-type  models,  while  for  other  classical  yield  functions, 
irrespective  of  the  procedure  followed,  the  constitutive  functions  a  and  i  cannot  be 
determined  or  the  resulting  kinematic  constraints  are  physically  unacceptable. 

We  shall  suppose  in  the  following  that  a,  0  depend  on  D  through  its  two  first 
invariants 

d^ir(D),  e  =  |XJ'|f  (9) 

i.c.  a  —  a((Le,p),{1  —  /?{d, c,p),  while  the  flow-no  How  function  F(a.p)  depends 
on  cr  through  the  invariants 

(10) 


(ID 

2,3  Per zyna- type  fluids 

A  general  procedure  for  extending  any  inviscid  elastic/  plastic  model  such  as  to 
account  for  rate  effects  was  proposed  by  Pcrzyna.  According  to  Porzyna’s  theory 
(1966),  viscous  plastic  flow  follows  a  time-dependent  How  rule  and  occ  urs  only  if  the 
stress  exceeds  the  yield  limit. 

Let  <t>  be  a  scalar  function  such  that  0(x)  =  0  for  x  <  0  while  fur  x  >  0,  we 
have  0(x)  >  0  and  <j>  is  monotonically  increasing.  Denoting  by  r;  >  0  the  viscosity 
coefficient,  we  consider  the  following  representation  of  the  rate  of  deformation 

°  =  +  B(<T,p)tr'  +  C(tr,p)(Til  ,  (12) 

where  F  is  the  yield  function  and  A,  B  and  C  are  scalar  functions  which  depend 
on  the  stress  invariants  and  on  the  density  p  wliich  plays  the  rule  of  hardening 
parameter  (sec  last  section  for  a  discussion  concerning  other  choices  for  the  internal 
variables). 

An  important  class  of  materials  can  be  derived  by  assuming  the  existence  of  a 
viscoplastic  potential,  i.e.  the  existence  of  a  function  G(a.p)  such  that 

A(tr,p)I  +  B(cr ,  p)cr'  +  C(crfp)tr2  =  Is)  (13) 

ua 


p  =  ^-,  ^  y  2 l<T 

i.c.  F  —  F(p.  p).  From  (65)  follows  that  for  D  ^  0: 

Q  =  J 1 0(d-  e,  p)e,  p  =  a(d,  e,  p). 
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Generally,  it  is  also  assumed  that  flow  is  associated,  Lc.  that  G  coincides  with 
the  yield  function  F .  However,  for  granular  or  porous  materials,  the  experimental 
evidence  suggests  that,  the  hypothesis  of  flow  normality  docs  not  apply,  and  the  more 
general  representation  (12)  of  the  rate  of  deformation  ought  to  be  used. 

We  seek  to  obtain  a  model  of  type  (65)  from  (12),  i.e.  to  deduce  explicit  ex¬ 
pressions  for  the  constitutive  functions  a  and  0  as  well  as  the  kinematic  restrictions 
resulting  from  enforcing  the  compatibility  conditions  (5)  to  (8).  Hence,  in  (12)  we 
neglect  the  influence  of  the  third  stress  invariant  and  take  C  =  t)  and  further  assume 
that  the  constitutive  functions  ,4  mid  D  as  well  as  the  viscosity  coefficient  tj  depend 
on  the  actual  density,  i.e. 

D  -  TTT <£(F(p,<7,p)}  \A(p,q,p)I+  B(p,q,p)(r'].  (14) 


Using  15and(14),  we  obtain  for  F(p,q,p)  >  0  the  following  nonlinear  system  in 
the  unknowns  p  and  q 


-<t>{F(p,q,p))A{p,q)  =  d 


2  nip) 

T-^\j^<i>(F'(p,q,p))qD(p,q)  =  e. 


2  ri{p) 


(15) 


Let  us  denote  by  (■,  •,  p)  :  R  x  R+  — >  1  x  R+,  the  function  defined  as  (e,d)  = 
SF{p,q,p)  through  (15).  Let  J$F  be  the  jacobian  matrix  of  SF.  Its  determinant  is 
equal  to  3\/2/3J/(16tj(/>)2)  with  J  given  by 

=  (<t>'AdpF+<t>dpA)(q<t>'Bdl,F+q<t>d<lB+<t>B)-{<i>'Ad,1F+<t>dllA)(<i<j>'BdpF+q<t>dpB) 

(16) 

For  a  particular  choice  of  the  constitutive  functions  A,  B  and  F,  the  set  of 
admissible  rate  of  deformation  D  is  given  by  the  image  of  the  operator  SF.  i.e. 


C  =  {D  €  R^*3;  3(p,  <7)  with  F(p,q,p)  >  0,  such  that  (d,e)  ~  SF(p,q,p)},  (17) 


while  the  invertibility  of  Sp  depends  on  the  rank  of  JsF  (given  by  J)  mid  dim(C). 

If  Sp  is  invertible,  i.e.  for  all  (d,  e)  6  C  we  can  solve  the  system  (15)  and  there 
exist  P  and  Q  such  that 


then 


at{d,e,p)  =  P{d,e,p), 


P(d,e,p) 

Q(d.  e,  p), 

(18) 

PL 

to 

1! 

ft  \ 

to 

(19) 

Next,  these  expressions  for  a  and  /?  are  substituted  in  the  compatibility  restrictions 
(5)-(8). 


9 


If  the  compatibility  conditions  arc  not  satisfied,  it  means  that  it  is  not  possible 
to  construct  a  fluid-type  model  (65)  using  the  regularization  method  of  Perzyna  (Lo* 
from  the  model  (14)), 

It  is  to  be  noted  that  if  the  hypothesis  of  associated  flow  rule  applies  and  0(;r)  — 
([x]+/Fo)n,  where  [x]+  —  (|x|  +  x)/2  and  Fq  and  n  arc  constants,  i,c. 


then  (15)  is: 


with  H(py  q,p)  —  (F(p,(p  p)/ F0)n+] /(n  +  1).  In  this  case  the  Jacobian  (16)  has  a 
simple  form 


3*H,  ,82H 


d2u 


J(p,q,p)  =  2F$[-j^[p,q,p)-^{p,q,p)  -  (^^(p,g,p))2|,  for  F{p,qrp)  >  0. 

{21) 


It  is  worth  noting  that  if  F  is  affine  in  p  and  q  (e.g.  Miscs  or  Di  ucker-Pragcr  yield 
functions)  then  J  is  always  vanishing  and  dim(C)  <  6,  i.c.  the  flow  has  to  satisfy 
certain  kinematic  restrictions.  However,  the  resulting  kinematic  constraints  have 
a  clem'  physical  meaning  only  for  certain  yield  conditions.  In  the  next  subsection, 
we  show  that  in  the  case  of  the  von  Miscs  yield  condition,  the  kinematic  constrainl 
that  ought  to  be  fulfilled  is  the  incompressibility  condition.  The  use  of  a  Drucker- 
Prager  yield  condition  would  lead  to  kinematic  restrictions  that  do  not  have  physical 
meaning  (see  section  5  for  a  detailed  discussion), 

2,3,1  M  ises  fluid 

Consider  the  Miscs  yield  condition  F{p.q)  —  \j\q  ft,  where  k  is  the  flow  stress  of 
the  material  in  simple  tension  or  compression.  For  associated  flow,  the  system  (20) 
reduces  to 


The  Jacobian  J  is  zero  (see  (21))  and  the  image  of  the  operator  Sjr  is: 

C  =  {D  €  K|x3  ;  trace(D)  =  0}, 


(22) 
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i.c.  the  flow  ought  to  be  isochoric.  Since  di.m{C)  —  5,  from  (20),  we  cannot  determine 
In  other  words,  p  cannot  be  determined  from  a  knowledge  of  the  deformation 
because  of  the  restriction  of  material  incompressibility.  From  (20)  follows  that  q  = 

Q(e,d)  =  ^|[k  +  Ff)(2?/(p)e/F1o}1/'!]  and  (65)  becomes: 


div(v)  —  0,  i 


er  =  al  4- 


*  +  F0(2r!{p)\D{v)\/F0)'/” 
\D(v)\ 


D(v)  if  |£>(v)|  *  0, 


W\  <  k  if|U(«)|  =  o, 

(23) 

with  n  undetermined.  Equation  (23)  can  also  be  viewed  as  a  constitutive  equation  for 
an  incompressible  viscous  fluid  with  viscosity  equal  to  («+Fo(2Tj(p)|jD(v)|/F0)l/T,)/!D(n)j. 

For  n  —  l  and  k  constant,  (23)  reduces  to  the  three-dimensional  form  of  the 
incompressible  Bingham  model  due  to  Hohcnemser  and  Prager  (see  for  example, 
Prager  and  Hodge  (1968)): 


div(v)  =  0, 


er  =  al  +  (2  i]{p)  +  |^^-)P(u) 
|<t'|  <  K 


i£\D(v)\^0t 

(24) 

if  \D(v)\=Qf 


The  Bingham  constitutive  equation  (24)  with  flow  stress  k  constant  is  the  most 
commonly  used  was  used  extensively  to  model  the  rate  and  time  effects  on  the 
deformation  of  metallic  materials  (  c.g.  Batra  and  Lin,  1987).  Bingham  model 
is  also  the  most  commonly  used  model  to  describe  fluids  with  yield  stress  such  as 
slurries  (  e.g.  aluminum  slurries)  and  waxy  crude  oils  (sec  for  example,  Alexandrou 
ct  ah,  2003;  Vi  nay  et  ah,  2005). 

However,  a  realistic  description  of  the  behavior  of  geological  materials  cannot 
neglect  the  dependence  of  yielding  and  subsequent  flow  on  density.  To  account  for 
this  influence,  Cristescu  et  al  (2002)  have  used  for  the  description  of  slow  motion  of 
natural  slopes  a  Bingham  model  (24)  with  yield  limit  dependent  on  density  k  =  s(p)< 
The  main  limitation  of  such  a  model  is  that  the  deformation  is  isochoric,  which  is  a 
direct  consequence  of  the  form  of  the  yield  condition  adopted. 


2,3-2  Mises  fluid,  revisited 

Implicit  in  the  models  of  type  (23)  is  the  assumption  of  incompressibility.  Also,  such 
models  cannot  account  for  changes  in  density  during  loading.  However,  experimental 
studies  of  the  dynamic  behavior  of  cementitious  materials  have  indicated  a  strong 
dependence  of  yielding  and  subsequent  flow  on  density  (see  for  instance  experiment  al 
data  on  cementitious  materials  reported  in  Malvern  et  al.T  1990;  Ross  ct  al.,1996; 
Schmidt  ct  al.T  2001  as  well  as  Hugoniot  data  reported  by  Osborn,  1981). 

To  describe  the  high-pressure  and  high  strain  rate  behavior  of  such  materials, 
we  propose  an  extension  of  the  Bingham  model.  We  assume  that  yielding  under 
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deviatoric  conditions  is  of  Von  Mises  type  with  yield  limit  k  depending  on  the 
current  density,  he.  k  —  n(p).  For  hydrostatic  conditions,  the  response  is  described 
by  a  pressure-density  relationship  of  the  form:  p  =  pc(p).  The  specific  expression  of 
the  function  pr[p)  is  obtained  from  quasi-static  and  dynamic  laboratory  data  as  well 
as  Hugomot  data.  Since  during  unloading  the  reversible  decrease  of  volume  is  very 
small,  it  can  be  neglected  and  therefore  the  unloading  process  is  rigid  (the  density 
acquired  at  the  end  of  the  loading  process  is  preserved  during  unloading-},  Thus,  in 
order  to  capture  both  the  effects  of  compaction  on  shear  flow  and  changes  in  density 
during  hydrostatic  loading,  we  consider  a  yield  condition  of  the  form: 


F(p,Q,p)=  \j -«(/»)]+  +  [p  +  pr(p)l_.  (25) 

where  [x]+  =  (|x|  +  x)/2,  [xj_  =  (\x\  —  x){2.  Next,  the  model  is  constructed  using 
the  procedure  outlined  in  section  3.1.  Indeed,  the  system  (20)  becomes 


^—{F(p,q,f))/Fa)n-llyj^q  -  «(p)]  + 


'Mp) 


ft 


Note  that  it  can  be  inverted  for  d  <  0  or 


C  —  {D  :  trace(D)  <  0}. 


These  kinemat  ic  restrictions  express  that  the  flow  is  compressible  and  the  unloading 
is  rigid.  Furthermore, 

,  -  PM ,)  -  ~pM + fot  ,  <  o 

V  a  +  e£ 

,  ,  F^{p)y[&T7*lF^I'\  r 

q  =  Q(d,  e,  p)  =  y  -[ K[p)  +  c - ~-j==- - 1,  for  e  /  0. 

Further  substitution  into  (65)  leads  to  the  following  constitutive  equation  of  a  com¬ 
pressible  rigid-viscoplastic  fluid: 


€T  — 


F0(2V(Pydiv(v)i  +  |D'|VF0)>A‘  «(/;) 

- r 


.  k'|  <  k(p) 
div(v)  <  0 


s/dMyFTW? 


IjD'I 


D'  if  \D'\  =£  0, 
if  \D'\  =  0, 


(2G) 


P  >  -Pr(p) 


+  |£>'|3 
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if  (/iu(-u)  —  0, 
(27) 
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Figure*  1:  The  dependence  of  pc  (left)  and  of  the  yield  limit  k  (right)  with  respect 
to  the  volumetric  strain  p  =  p/pa  —  1  for  concrete  (data  after  [30]). 


This  model  is  an  extension  of  Bingham  type  model  (23)  that  captures  the  experi¬ 
mentally  observed  density  changes  under  severe  loadings.  If  we  set  n  =  1,  then  we 
obtain: 


a  — 


2tj(p)  + 


k(p) 


\D'\\ 


D'  if  \D'\  ^  0, 


(28) 


W\  <  k(p) 


if  |D'|  =  0, 


{trac.e(c r)/3  =  -pc(p)  +  2r}(p)div(v)  if  div(v)  <  0 

(29) 

froce(er)/3  >  —pc(p)  if  div(v)  —  0. 

In  what  follows,  the  constitutive  equation  (28)  will  be  referred  to  as  the  ’’revisited 
Bingham  fluid”. 

In  many  solid-type  models,  and  k  arc  expressed  as  functions  of  the  compaction 
level  (or  volumetric  strain)  p  —  p/pq—  1  instead  of  the  density  p  (here  po  is  the  density 
prior  to  loading).  As  an  example,  in  Figure  1  is  plotted  a  typical  dependence  of  pc 
and  of  the  yield  limit  k  versus  the  compaction  level  p  for  a  concrete  material  (data 
mid  explicit  expressions  of  these  functions  reported  by  Holmquist  ct  al.,  1993),  The 
explicit  expressions  of  these  functions  arc  given  in  section  5.  Furthermore,  the 
revisited  Bingham  fluid  model  can  be  recast  in  the  form  that  is  generally  used  in 
hydrocode  calculations,  by  eliminating  p  between  k  =  n{p)  and  pc  =  pc(p)-  This 
type  of  dependence  is  also  used  in  soil  mechanics  to  eliminate  the  density  as  an 
unknown  function.  For  the  same  concrete  material,  the  obtained  relation  between 
the  shear  flow  stress  and  pressure  is  represented  in  Figure  2. 
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Figure  2:  The  dependence  of  of  the  yield  limit  k  with  respect  to  the  pressure  p  for 
undamaged  concrete  (data  after  Holmquist  et  ah,  1993). 


2.3.3  Modified  Cam-Clay  fluid  of  Perzyna’s  type 


As  already  mentioned,  there  is  a  critical  need  for  developing  constitutive  models 
that  could  capture  the  complex  rheological  behavior  of  slurries  and  some  oils.  The 
fundamental  question  which  arises  concerns  the  (low  pattern,  particularly  the  shape 
and  distribution  of  yielded /unyielded  regions. 

The  Modified  Cam-Clay  yield  function  developed  by  ftoscoe  and  Garland,  1908 
captures  very  well  commonly  observed  properties  of  soils  such  as  compact  ion/dihitaney 
behavior  and  the  tendency  to  eventually  reach  a  state  in  which  the  volume  becomes 
constant.  It’s  adequacy  has  been  tested  against  numerous  experiments,  the  consen¬ 
sus  being  that  is  most  suitable  for  lightly  overconsolidated  particulate  materials  (for 
a  review  on  the  subject,  see  dc  Burst  and  Croon,  2000). 

In  view  of  applications  to  the  description  of  the  flow  of  slurries,  we  further 
examine  whether  a  compressible  viscoplastic  fluid  model  compatible  to  a  Modified 
Cam-Clay  yield  condition  can  be  constructed  using  the  viseoplastie  regularization 
method. 

The  Modified  Cam-clay  yield  com  lit  ion  is 

F{p,  q ,  p)  =  \/q2jM2  +  (p  +  pc{p))2  -  pe(p), 


wlicrc  A/  is  a  constant.  If  we  assume  associated  flow,  the  system  (20)  becomes 


Fo 


{ F(p,q,p)/F0)n 


V  +  Pc{p) 


F0 


.  2r)(p) 


2v(p) 
\jl{F(p,q,p)/F0)n 


yV/A/2  +  (p  +  pr(p))2 
_ q/M2 _ 

\J  q1  j  M2  +  (p  +  Pr{p))2 


=  d 


=  e. 
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In  contrast  to  the  above  example  (Mises  yield  condition),  the  system  can  be  solved 
for  all  (d,  e)  £  R  x  R+,  i.c.  there  arc  no  kinematic  restrictions  C  =  After  some 

algebra,  wo  obtain: 


P(d,etp)  =  -pt(p)  +  d 


P*{p)  +  F0(  Mph/^  +  2MV/3/fB)‘^ 

y/tP  +  2A/2ea/3 


.  e.  P)  =  \/|M2fi 


pc(p)  +  F0(2r}{p)y/tP~+2M^yi/F0)ifn 


Q{d,~,ri  —  m  -  , - 

l/  "  +  2A/2e2/3 

Substitution  in  (65)  leads  to  the  following  constitutive  equation: 


cr  = 


~Pr[p)  +  div(v) 


Pc{p)  +  FQ(2T}{p)^/div(v)2  +  2M2|D'(u)|2/3/F0))/n 
\AMo)a  +  2A/2|D'(v)|2/3 


/  + 


2A/aPe(p)  +  fo(2q(p)>/*t;(i>)a  +  2Ma|C'(t>)|2/3/F0)1'"  ,f  |n,#oj  ^  n 

3  v/^(u)2  +  2M2|D'(u)|V3 

F(p,?,p)  <  0  if  |D(w)|  =  0 
(30) 

The  constitutive  equation  (30)  describes  a  compressible  fluid  with  flow  behavior 
dependent  on  the  first  stress  invariants  and  strain  rate.  If  we  set  n  =  1,  we  obtain 
the  following  simplified  model 


a  = 


-pc{p)  +  div(v)  (  2p(p)  + 

2  A/2 


Pc(p) 


2 n(p)  + 


y/div[v)2  +  2A/2|Z3 

Pc(p) 


1+ 


_ _ T _  D'(v)  if  |jD(®)|  ?  0, 

\/div(v) 2  +  2A/2|.D  (v)|2/3 

F{p,q,p)  <  o  if  I £»(w)l  =  o. 

(31) 


2.4  Superposition  method 

In  this  section,  we  introduce  a  different  procedure  for  determination  of  the  explicit 
expressions  of  the  functions  a(d,  e,p)  and  0(d,e,p)  involved  in  the  general  constitu¬ 
tive  equation  (65),  The  main  assumption  is  that  the  state  of  stress  in  the  material, 
cr.  can  be  represented  as  the  sum  of  a  viscous  contribution  <rh  and  a  contribution 
S.  related  to  plasticity  effects,  i.c.: 

<r  =  aF  +  S.  (32) 

The  viscous  part  of  the  stress  is  taken  as  for  a  classical  viscous  fluid, 

aF  =  fv(D,p).  (33) 
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It  should  be  noted  that  at  difference  with  the  stress-rate  of  deformation  relation  for 
a  rigid  viscoplastie  fluid,  the  constitutive  function  /l  is  continuous  and  vanishes  for 
D  =  0,  lc. 

fv(D,p)^Q  =  fv(Q,p),  for  D  — *  0,  (34) 

For  the  sake  of  simplicity  we  shall  consider  only  first  order  fluids 

fV  (D.p)  =  \{D,p)tr{D)I  +  2t ]{D,p)D,  (35) 

where  T)  and  A  arc  viscosity  coefficients  which  may  depend  on  d,e  mid  p*  i.e.  A  = 

A (rf,e, /?),//  —  r)(d,  e.p),  For  A  and  r/  independent  of  d  and  e%  i.e.  A  =  A(/j)t7j  =  q{p), 
(33)  reduces  to  the  compressible  Navicr-Stokcs  model 

fv(D,p)  =  \(p)dI  +  2r,(p)D ,  (30) 

but  other  choices  can  also  be  considered  in  (35),  For  example,  for  q  —  B(p)  arg  sinh(j4e)/ 
with  A,B  >  0,  we  recover  the  P  rand  tl*Ey  ring  type  model  and  for  //  —  fi(p)ca,  with 
a  =  1/m  —  1  >  —1,  we  obtain  the  Norton  model.  Models  of  isotropic  fluids  have 
been  used  to  describe  the  slow  motions  of  soils  and  glaciers  on  natural  slopes  (sec 
Vullict,  2000  for  an  overview  on  the  subject). 

We  assume  that  there  is  flow  only  if  a  yield  condition  expressed  in  terms  of  S 
and  depending  on  is  satisfied  (for  a  discussion  concerning  other  choices  of  the 
hardening  or  damage  parameters  see  last  section).  We  further  neglect  second  order 
effects  in  S  as  well  the  influence  of  the  third  invariant  of  S.  Thus,  the  yield  condition 
is  expressed  as  F(s ,  r  p)  —  0,  where  s  —  lr(S)/3  and  r  —  y/3/2  |S'|  arc  the  first  two 
invariants  of  S.  The  deformation  rate  D  is  then  given  by 

D  -  A  [M(a9r,p)I  +  N(s,r,p)Sf  ] ,  (37) 

with  F(s,  t\p)  <  0  ,  A  >  0  and  \F(s,  i\p)  =  0,  while  M  and  N  are  scalar  functions 
of  their  arguments.  From  (37)  follows  that  for  D  ^  th  we  obtain  the  Following 
nonlinear  system  in  the  unknowns  s,  r  and  A 

f  3A  M(a,r1p)=  rf 

<  ArAT(s,7'1p)  =  e  (38) 

*  F(s,r,/>)«  0. 

Let  us  denote  by  TZf{%  *,  p)  :  R  x  K+  x  R+  R  x  1+  the  function  defined  by 
(38),  i.e.  (c,rf)  ~  Avp).  The  set  of  admissible  rates  of  deformation  D  is 

given  by  the  image  of  the  operator  i.e, 

C  ~  {D  £  R|x3;3(s,r,  A)  with  F(s1r,p)  =  0,  such  that  (d,e)  =  lZjr[s,  rT  A,p)}. 

(39) 
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The  invertibility  of  7 tf  depends  on  the  rank  of  the  Jacobian  matrix  of  Ti?  and 
dim(C).  If  (38)  is  invertible,  i.e.  for  all  (d,  e)  €  C  we  can  solve  the  system  (38)  and 
there  exist  5,  R  and  A  such  that 


{s=  S[d,e,p) 

r  =  R(d,e,p)  (40) 

A  =  A  (d,e,p), 

by  combining  (40),  (32)  and  (35,  we  can  obtain  the  explicit  expressions  for  the 
constitutive  functions  a  and  0  of  (65): 

a(d,  e,  p)  =  (A  (d,  e,  p)  +  Cl  +  S{d,e.p), 

<  r-'  ,  (■'ID 

a/ 1  v  „  s,  /2  R(d,e,p) 

P{d,  e,  p)  =  2r](d,  e,  p)  +  W  - - - - •• 

As  in  the  ease  when  we  make  use  of  the  Perzyna  viscoplastic  regularization 
method,  we  substitute  the  expressions  (41)  in  the  compatibility  conditions  (5)-(8). 
If  one  of  these  conditions  is  violated,  we  deduce  that  we  cannot  construct  a  fluid- type 
model  (65)  using  the  superposition  method  (32)-(37). 

In  the  following,  we  construct  fluid  models  using  the  superposition  approach  and 
compare  their  response  with  the  viscoplastic  fluid  models  presented  in  section  3. 
In  some  examples  we  shall  assume  an  associated  flow  rule,  i.e.  that  there  exists  a 
function  F(s,r,  p)  such  that 

M(S,p)I  +  N(S,p)S'=?£^-. 

In  such  cases,  the  scalar  scalar  functions  M,  N  arc  expressed  as: 

i  a  rp  o  i  j  p 

M(s,r,p)  =  ~—(s,r,p),  N(s,r,p)  =  —  ~^(s,r,p).  (42) 


2.4.1  Bingham  Fluid 

Wc  first  consider  the  superposition  between  a  viscous  fluid  and  a  plastic  solid  obeying 
Miscs  yield  condition,  i.e.  F(s,r,p)  =  r  -  n{p).  If  we  assume  associated  flow 
(42),  then  (38)  becomes 


This  means  that  the  flow  is  incompressible,  i.e.  the  kinematic  restriction  coincides 
with  (22).  Moreover,  s  and  a  arc  undetermined  while  R(d,e,p)  —  ^3/2 n{p). 
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If  we  choose  the  compressible  Navier-Stokes  model  (3G)  for  the  viscous  part  then 

j 3{d*ctp)  —  2 //  +  ”,  and  we  obtain  the  constitutive  relation  (24)  of  a  Bingham  fluid 
e 

with  yield  limit  dependent  on  density* 

Let  us  note  that  if  the  flow-no  flow  condition  is  Von  Miscs  yield  function,  using 
either  method  we  obtain  the  classic  Bingham  fluid  model  This  is  not  the  case  rot- 
other  yield  conditions. 


2,4,2  Bingham  fluid,  revisited 

Next,  we  consider  the  superposition  between  a  viscous  fluid  and  a  plastic  solid  that 
obeys  the  yield  condition 

F(s,r,p)  =  \y/2jlr  -  «(/>)]+  +  [s  +  pR(p)]_, 

where  [a‘]  +  =  (\x\  +  ar)/2,  [x]_  —  (\x\  —  x)/2,  and  associated  flow.  The  system  (38) 
becomes 

d  =  0,  e  =  A  if  s  >  -pc(p),  r  =  (p), 

d  =  —A,  e  =  0  if  s  =  —  pe{p),  r  <  \/3/2 s(p). 

If  s  ~  —pc{p)  and  r  —  y^3/2 *(p),  wc  have  to  consider  the  normal  cone  to  the 
plasticity  surface  and  M  and  N  arc  set  valued  functions,  i.e, 

d  <  0,  e  >  0,  if  s  —  —  Pc[p)*  t  —  y^3/2ft(p). 


The  kinematic  restriction  is  C  =  {Z?  ;  lr(D)  <  0}  ,  which  implies  that  the  flow 
is  compressible  and  that  unloading  is  rigid.  Furthermore,  for  |jD|  ^  0  the  system 
(38) can  be  inverted  and  wc  obtain: 


S(d:  e,  p)  =  ~pr(p)7  R{d,e,p)  < 


A (d.  e,  p)  =  —d  if  d  <  0,  e  =  0, 


S(d.e.p)  >  -Pc{p),  R{d 


•e'rt  =  vf 


«(p),  A{rf,  t,p)  =  e  if  d  =  O.e  >  0. 


Substituting  into  (41)  and  choosing  the  compressible  Navier-S  tokos  model  (36)  for 
the  viscous  part,  wc  obtain  the  following  constitutive  equation: 


1  »'= 

hw+Si] 

l  k'l  < 

K(p) 

D‘  if|D'[^0 
if  |£?'|  =  0, 


(43) 


div{v)  <  0, 


trace(tr )/3  -  -pc(p)  +  [A(p)  +  2ri(p) /3\div{v) 
trace{c r)/3  >  —pr(p) 


if  dh'(v)  <  0 

if  div(v)  =  0, 
(44) 
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Let  ils  note  that  this  model  is  different  from  the  constitutive  equation  (28)-{29) 
obtained  using  Perzyna’s  approach.  The  difference  stands  in  the  viscous  terms  of 
the  spherical  part  of  the  stress. 

Under  hydrostatic  conditions,  most  cementitious  materials  show  a  highly  non¬ 
linear  pressure- volumetric  strain  response,  the  reversible  decrease  in  volume  being 
very  small.  The  experimental  observations  also  suggest  that  in  the  high-pressure 
regime,  a  very  large  increase  in  pressure  is  necessary  in  order  to  produce  even  a  very 
small  change  in  density.  Thus,  the  hypothesis  of  a  ” locking  medium1*  can  be  made, 
he.  the  density  cannot  exceed  a  critical  value.  This  critical  density  p*,  called  locking 
density,  corresponds  to  a  state  in  the  material  when  all  the  pores  and  micro-cracks 
arc  closed.  The  pressure  level  at  which  this  density  is  first  reached,  called  locking 
pressure,  is  denoted  by  p *  =  pc(p*).  Hence,  for  such  materials  a  pressure-density 
relationship  of  the  following  form  could  be  considered: 


P  =  Pc(p),  it  p<  p* 

P  >p\  if  P  =  P‘- 


(45) 


2.4*3  Superposed  Cam-Clay  fluid 

Next,  we  consider  the  superposition  between  a  viscous  fluid  and  a  plastic  solid  that 
obeys  the  Modified  Cam-clay  yield  condition  F(s,r, /?)  =  y/r2/M2  +  ( s  +  pCr[p))2  — 
Pc[p)  and  an  associated  flow  rule.  The  system  (38)  becomes 


d  =  A 


S  +  Pc[p) 
Pc(p) 


-*  i/I 


JpjM'  jj; +  <•  +  *<*))’- rfM 


This  algebraic  system  can  be  solved  for  all  [d,  e)  €  lx  R+,  (i,e.  there  are  no 
kinematic  conditions  C  =  R|*3)  and  after  some  algebra  we  get: 


S{d,e,p)  =  -Pc(p)  +  d 


Pcip) 

^/d2  +  2M2e2/3’ 


R(d,e,p) 


f-M\  ,  ^ 

V  3  y/cp  +  2A/2e2/3 


Finally,  if  the  compressible  Navier-Slokes  model  (36)  for  the  viscous  part  is  chosen, 
we  deduce  the  constitutive  equation 


<7  = 


(»)  ^ 


-Pc(p)  +  div(v)  I  A(p)  +  + 


Pc{p) 


3  ^dn-'(v)2  +  2M'2\D‘  [v)\2 /Z  J  J 

D‘{v) 


I + 


2tj (p)  +  -  M2  Pr^ 

3  s/dMvj*  +  2A/2|D»|2/3_ 


if  \D(v)\  ±  0, 


,  F(p,q*p)  <  o 


if  |D(u)|  =0. 
(46) 

Let  us  note  that  the  above  constitutive  equation  for  the  Cam-clay  fluid,  obtained 
using  the  superposition  method,  is  slightly  different  from  (31)  which  was  constructed 


19 


using  Pcrzyna’s  approach  and  a  static  yield  condition  of  the  same  functional  form* 
In  both  models,  the  part  of  the  stress  which  is  rate  independent  is  the  same,  and 
only  the  viscous  contribution  is  different. 


2*4*4  CAP  fluid 


Cap  models  are  among  the  most  widely  used  models  for  the  description  of  the 
mechanical  behavior  of  soils*  We  examine  herein  whether  a  cap  yield  condition  can 
be  used  to  construct  a  viscoplastic  fluid-type  model  using  the  superposition  method. 
Thus*  we  superpose  a  viscous  fluid  on  a  plastic  solid  satisfying  a  cap  plasticity 
yield  condition  with  a  modified  Drucker-Prager  envelope  (sec  Katona,1984)1  i*e* 
we  consider  that  the  yield  surface  is  divided  into  three  regions  along  the  s-axis  (s 
being  the  first  invariant  of  the  plastic  component  of  the  stress  tensor):  the  failure 
surface  region  3s  €  (—L(p)>T);  the  cap  surface  region  3s  €  (-X{p),-L(p))  and 
the  tension  cut-off  region  (3s  >  7\  We  consider  the  failure  surface  to  be  of  a 
modified  Drucker-Prager  form  (see  Katona,  1984),  i*e*  rfyf 3  +  Oexp(3Bs)  —  ,4, 
where  A  and  C  are  constants  (A  >  (7).  It  constitutes  a  boundary  along  width  the 
cap  surface  can  move  (harden  or  soften)*  The  cap  surface  is  a  hardening  surface 
in  the  shape  of  an  ellipse  when  plotted  in  the  space  (r,  —  s)  (see  Figure  3.  It  is 


defined  as  j- 


73 


jX(p)  L(/i))Mft.^(p)) 

w 


The  variable  cap  parameters  L  and  X 


arc  positions  on  the  s  axis  which  locate  the  current  cap  surface,  while  the  parameter 
R  defines  the  ratio  of  the  principal  ellipse  radii.  The  tension-cut  toff  criterion  is 
triggered  when  s  >  1\  whereas  T  is  a  material  constant  representing  the  threshold 
of  hydrostatic  tensile  stress  at  which  abrupt  stress  release  occurs  due  to  tension 
damage*  Thus,  we  consider 


3,s  -  T, 


if  3,s  >  T 


F(s,r,p)  =  i 


rj  \/3  +  Ocxp(3Z?s}  -  A, 


To 


r2/3  - 


(X(p)  z  uytf  -  (3s  -  /.(;,)) 


if  -  L(p)  <  3s  <  T 
if  —  X(p)  <  3s  <  —£(/;), 


(47) 

and  an  associated  flow  rule  (see  a  depiction  of  the  yield  surface  on  Figure  3  (left  ). 
Then  the  system  (38)  becomes 


d  =  3A  e  =  0 

d  -  3Af?Ccxp(3Bs),  e  =  A/ ^2 

d  =  ~ Lip)l  e  =  \ZfNkr 


if  s  =  T/3,r  e  (0,rr), 

if  3s  €  {—L(p),T),  r=s/3{A  +  C  cxp{3 Bs) 
if  33€  {-X(p),-L(p)),F{s,r>p)=0 


(48) 

if  the  stress  point  (s,  r)  is  on  the  smooth  part  of  the  plastic  surface  F(s,i\p)  —  0. 
In  the  two  singular  corner  regions  3s  =  1\  r  =  /?■/■  and  3s  =  -L(p),r  =  Rl(p)  (see 


2t) 


figure  3  left)  the  gradient,  of  F  is  set  valued  and  we  have 


f  3ABC cxp(BT)  <  d  <  3A,  0  <  e  <  A/\/2  if  3s  =  T,r  =  RT, 

{  0<d<3ABCcxp(-BL{p)),  A/y/2  <  e  <  yfeh-j-Rdp),  if  3s  =  -L(p), r  =  RL(p) 

(49) 

The  above  algebraic  system  can  be  solved  for  all  (d,e)  e  R  x  Et+,  (i.c.  there  are 
no  kinematic  conditions  C  —  R|*3)  and  after  some  algebra,  we  get: 


r  S{d,e,p)  =  Tf 3,  R(d,e,p)  =  RT 
S{d,e,p)  =  gL[ln(d/e)  -  ln(3v^BC), 
R(d,  e,  p)  =  \/3[A  -  ^j^dj e] 


< 


S{dye,p)  =  ~L(p)/ 3,  R(d,e,p)  =  RL(p ), 


R(d ,  e ,  p) 


JML 

»o  s/rV3+R^JbA 


if  d/e  >  3\/2BC  cxp(BT), 

if  3\/2BC  cxp(BT)  >  d/e.  >  3\/2BCcxp{— BL(p)), 
if  3s/2BCcxp(-BL{p))  >  d/e.  > 

if  _  flJ*!AsL  >  d/e 
V  3  R%Rai)  ~  a/e- 

(50) 


Finally,  we  arrive  at  the  following  cap  fluid  type  model 


<T  = 


div( 


«)  ( 


2,!W  +  V  3 - iD’M - 


+  S(diy(w).  |D'(«)|,p)  r+ 
0'(«) 


,  B(p,(?.p)  <  0 


if  |z>(*)|^o,  (51) 

if  |i?(w)|«0, 


if  the  compressible  Navicr-Stokcs  model  (36)  for  the  viscous  pai't  is  chosen. 


2.5  It  does  not  always  work 

As  mentioned  previously,  it  is  not  always  possible  to  obtain  consistent  viscoplastic 
fluid  formulations  starting  from  any  two- invariant  yield  functions.  A  first  example 
is  given  below. 

Let.  F(p,  q.  p)  —  q+p  tan{/3)  — «o (p)<  i.c.  Drucker-Prager  yield  condition  with  con¬ 
stant  angle  of  friction  /?  and  cohesion  k()  dependent  on  the  density.  We  take  0(.r)  = 
([:r]+/Fo)n  and  further  assume  non- associated  (low,  i.e.  A(p.q,p)  =  FndpG[p,  q)/3 
and  B(p.q)  =  3Ff}d„G{p,q)/{2q).  It  follows  that 


F  BO 

d=m{lF{p'q’p)h/F°rwM' 


=  4 


r 


If  the  potential  is  chosen  of  the  form  G(/M/)  =  9  +  ptan(^),  where  the  constant  ift  is 
the  dilation  angle,  then  from  the  above  equation,  we  deduce  that  d  —  ^/2/3  tan(^)e. 
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Figure  3:  Schematic  representation.  Loft:  the  cap  plasticity  condition  with  a  mod¬ 
ified  Druckcr-Prager  form  used  for  the  cap  fluid  (51).  Right:  the  cap  plasticity 
condition  with  a  standard  Druckcr-Prager  surface  which  violates  the  continuity  con¬ 
dition  (6),  The  shaded  areas  are  ” singular  corner  regions*' . 

This  means  that  the  set  of  admissible  deformation  rates  is: 

C  =  {D  €  K|*3  ;  tr{D)  =  tan(^)|D'|}.  (52) 

V 

This  kinematic  constraint  is  rather  artificial  and  severely  limits  the  flow.  Moreover, 
reinforcing  such  a  constraint  would  pose  insurmountable  numerical  problems. 

If  we  use  for  the  flow  potential  another  wcll-acceptcd  expression  (e.g.  see  [2|) 
G(p,q)  =  \J V2  +  Po  +  ptan(^),  where  Iq  and  are  constants,  we  obtain: 

d  = 

which  implies  the  following  kinematic  admissible  set 

C  =  {D  6  K3*3  ;  tr(D)  >  tan(^)|I?,|}.  (53) 

As  in  the  previous  case,  the  resulting  kinematic  constraints  restrict  drastically  the 
flow  and  arc  very  difficult  to  justify  or  handle  numerically. 

Concerning  the  Cap  fluid  model  (51),  it  is  worth  noting  that  even  though  the 
plasticity  condition  is  specified  in  terms  of  three  functions,  the  inverse  functions 
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S  and  R  arc  continuous  in  R  x  R+  \  {(0,0)}.  This  is  a  direct  consequence  of  the 
choice  of  an  exponential  dependence  on  the  first  invariant  instead  of  the  standard 
Drucker-Pragcr  expression  of  the  envelope.  Next,  we  show  that  if  the  standard 
D  rocker- Pragcr  envelope  is  used  then  the  response  functions  R,S  are  discontinuous, 
which  is  physically  unacceptable. 

Indeed,  if  in  the  cap  plasticity  condition  (47),  with  associated  flow  rule,  we  take 
(sec  [20])  the  standard  Druckcr- Pragcr  expression  for  the  envelope,  i.c.: 

F(s,  r,  p)  =  r/y/Z  +  ZDs  -A  if  -  L(p)  <  3s  <  T.  (54) 


(see  also  Figure  3  (right)  then,  after  some  algebra, we  find  that  there  are  no  kinematic 
conditions,  C  —  R^*3,  and  wc  obtain  the  following  expressions  of  S  and  R  instead 
of  (50): 


S{d,  e.,p)  =  Tj 3,  R[d,e,p)  =  RT 
S(d,e,p)  =  -L{p)/ 3,  R(d,e,p)  =  RL(p) 
S(d,e,r)  =  L(P)/ 3+no^l^lli-, 


if  d/e  >  3  y/2/B, 
if  Zy/2/D  >  d/e  > 


i  »/-(p) 


l(p)' 


if 


1  »J4p) 


>  d/e 


(55) 

Note  that  the  functions  S  and  R  have  a  discontinuity  on  the  line  d  —  3\/2e/D  of  the 
d  —  e  plane.  From  (41)  it  follows  that  the  constitutive  functions  a  and  have  the 
same  discontinuity  with  respect  to  D,  which  violates  the  compatibility  condition 
(6). 


2.6  Examples  of  the  models’  response 

Wc  illustrate  the  response  of  the  models  proposed  for  uniaxial  strain  compression 
conditions.  Suppose  that  in  the  reference  configuration  the  material  sample  occu¬ 
pies  the  cylinder  Co  =  (0.  £<})  x  u;,  In  the  deformed  configuration  at  time  /,  the 
material  sample  occupies  Ct  —  (0.  £(£))  x  uj  (see  Figure  4  top).  The  deformation 
is  supposed  to  be  homogeneous  x  —  L(i)X/Lv,  y  =  Y f  z  —  Z  and  the  velocity  in 
Euler  ian  coordinates  xiy1  z  is  vx  —  L(t)x/L(t)y  vy  —  0,  vz  =  0,  so  the  only  non-zero 
component  of  the  rate  of  deformation  tensor  is  Dxx  =  L(t)/L(t). 

If  we  denote  by  po,  the  density  at  t  =  0,  then  at  any  time  t  >  0,  the  density  p , 
the  compaction  factor  (volumetric  strain)  fi  =  p/p®  —  1,  and  the  invariants  of  D  arc 
given  by 


,  ,  L(t) 
Pit)  =  -r^P°’ 


Hit)  = 


La 

L(t) 


-  1, 


/m  - 

d(  )  ~  3  L(t)' 


p\kt)\ 

e(f)  =  VsW 


This  loading  process  corresponds  to  compression  of  the  sample  material  in  the  ;r 
direction  (see  Figure  4  top). 
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Figure  4;  The  evolution  of  the  compaction  factor  t  — *  p(t)  (left)  and  of  the  equivalent 
strain  rate  t  — e(Z)  (right)  for  the  second  loading  history. 


Let  us  consider  the  following  loading  history:  the  piston  which  at  t  —  0  has  the 
velocity  V^Q)  =  0,  uniformly  accelerates  in  the  time  interval  t  6  [(Mej  to  reach  a 
constant  velocity  value  V,  which  is  kept  constant  for  t  e  \toAu  +  71],  i.c.  K(£)  =  V 
and  then  uniformly  decelerates  to  reach  a  vanishing  velocity  at  /  =  T  +  2 4-  That 
means  that  the  velocity  of  the  piston  is  given  by  L(t)  —  —  iV/to  for  t  E  [0,lo], 
L(t)  —  — V  for  t  E  [foT  Iq  +  T]  and  L{t)  —  —  V  ( T  +  2£q  —  t) /t^. 

To  examine  the  capability  of  the  proposed  models  to  describe  high-strain  rate 
effects,  we  compare  the  response  of  each  material  for  three  different  velocities: 
V/Ld  =  Ws~]  ,V/Ld  —  1 00s"1  and  V/Lq  =  1000s"1,  We  choose  the  constants 
Iq  and  T  such  as  to  ensure  that  the  final  compaction  level  pjinai  —  0.06  reached  in 
each  experiment  is  the  same  (i,c.  we  take  tG  =  1/240$ fT  =  1/80 s  for  the  first  ex¬ 
periment,  =  l/2400stT  —  l/800s  for  the  second  and  —  1  /24000s,  T  =  1 /8000s 
for  the  third  one).  As  an  example,  the  evolution  of  the  volumetric  strain  /  — ►  /i(f) 
and  of  the  equivalent  strain  rate  t  — ►  e(£)  for  the  second  loading  history*  are  plotted 
in  Figure  4,  Note  that  in  the  Eulerian  description  the  strain  rate  is  not  constant  on 
t  €  [^o,io  +  T]  even  if  the  piston  has  a  constant  velocity. 

2*6,1  Bingham  fluid,  revisited 

First,  we  consider  the  revisited  Bingham  fluid  model  (G8}~(69).  The  specific  expres¬ 
sions  of  the  constitutive  functions  k(p)  and  pdp),  taken  from  Holmquist  et  al.  (1993), 
correspond  to  a  concrete  material  (see  also  Figure  1  for  plot  of  these  functions).  The 
function  p?  is  expressed  in  terms  of  the  volumetric  strain  p  =  p/po  -  1  as:  pr{p)  — 

PP crush  / Pcrush  if  P  5:  Pcrush  i  Pr(p  )  “  Pcruch  T  {p  P crush  )  (PforA  Pcruch )  /  ( Ptock  Pcrush  ) 
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Figure  5:  Left:  the  stress  evolution  /,  — *■  (— p(t),ff(0)  of  the  revisited  Bingham 
model  (68)- (69)  for  three  loading  histories.  Right:  a  comparison  between  the  stress 
evolution  (  — ♦  (— p[t),q(t))  for  the  two  revisited  Bingham  (BR)  models,  (28)-(29) 
obtained  using  the  Perzyna  approach  and  (68)- (69)  obtained  using  superposition 
approach,  and  for  V/L0=  1000/s. 
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Figure  6:  Right:  the  stress  evolution  f  — >  (-p(/),  r/(0)  of  the  Cam-clay  fluid  (31) 
for  three  loading  histories.  The  solid  line  represents  the  final  plastic  surface,  i.c. 
F(p,q,  =  0.  Left:  the  very  beginning  of  the  stress  path.  The  solid  line  repre¬ 
sents  here  the  initial  plastic  surface,  i.c.  F(p,q,  0)  —  0. 


if  M  £  [Pcmsh,p.\  and  pc(p)  =  K\ {p- filock)/[l  +  W<**)]2  + 

-  plocjt)/(l  +  fkock)}A  if  P  >  p.-  The  specific  values  are  =  16  MPa, 

Pior jt— 800  MPa,  fi crush—  0.001,  ftiock  =  0-1,  //*  =  o-lll  and  A'j  =85000  MPa,  /v*=- 
171000  MPa,  A'j=208000  MPa.  The  function  K(fi)  depends  on  ft  through  pc(p)  (see 
Figure  2) 

n(p)  =  oQ\j2]Zmm{Smai,  A  +  B{pc(p)/trQ)N }  (56) 

with  <7o=48  MPa,  Smax  =  7,  A  =  0.79,  B  =1.60.  N  —  0.61.  The  viscosity  coefficients 
were  supposed  to  he  constants  A  =  0*002  MPas,  ??=  0.001  MPas. 

In  Figure  5  (left)  is  shown  the  evolution  of  the  stress  in  the  q  —  p  plane  for  the 
three  loading  rates  considered*  For  the  test  at  (Vr/£f(J=10/s)  and  (V/ A 0=100/*). 
the  stress  curves  arc  very  close  to  the  the  static  curve  q  —  \j3/2t\{p),  A  larger 
difference  with  respect  to  the  static  response  is  observed  in  the  third  experiment 
[V/L0=  W00/s). 

A  comparison  between  the  revisited  Bingham  (BR)  models  obtained  using  the 
Perzyna  approach  (28)-(29)  and  the  superposition  method  (68)- (69)  is  presented 
in  Figure  5  (right)  for  V/Lq=1000/s ,  The  rate  effects  as  well  as  the  difference  in 
response  between  the  two  models  is  clearly  observed* 
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2.6.2  Cam-clay  fluid 


Wc  apply  the  Cam-clay  fluid  (31)  obtained  using  Pcrzyna  approach  to  a. shale  (data 
after  Niandou,  1994),  The  expression  of  pc  is:  pr(p)  —  po(l+/*)1^  which  corresponds 
to  an  oily  shale  sample  (Tournemire,  France).  The  material  constants  arc:  M  = 
2,6,  po  =  2  MPa,  l— 0-02  and  the  viscosity  coefficient  is  supposed  to  be  constant  ij 
=  0.01  MPas* 

In  Figure  6  (right)  we  have  plotted  the  evolution  of  the  stress  in  the  q—p  plane  for 
the  three  loading  rates.  The  final  plastic  surface,  i*c.  F(p,<j,  (ifin)  =  0  is  represented 
by  a  solid  line.  For  the  first  two  loading  rates  V/L{)  —  10/s  and  V/L fl  =  100/s, 
the  response  is  very  close*  Some  difference  in  response  is  observed  for  an  order  of 
magnitude  increase  in  the  Loading  rate  ( V/Lq  =  1000/s).  Figure  6  (left)  shows  the 
beginning  of  the  loading  process.  In  this  figure,  the  solid  tine  represents  the  initial 
plastic  surface,  i.c.  F{p,q,  0)  —  0.  The  rate  effects  on  the  flow  are  clearly  observed. 
Although  all  the  stress  curves  start  from  the  initial  plastic  surface,  their  evolution  is 
different,  the  most  pronounced  difference  being  between  the  response  corresponding 
to  the  lowest  and  highest  loading  rates. 

2,7  On  the  choice  of  internal  parameters 

In  all  the  models  presented,  the  assumption  was  made  that  the  only  dissipative 
mechanism  is  plastic  flow  and  the  associated  internal  variable  is  either  the  density 
p  or  the  porosity/compaction  level  p.  Although  the  choice  of  the  density  as  hard¬ 
ening  parameter  is  physically  sound  and  very  suitable  for  an  Eulerian  description 
of  large  deformation  and  high  strain  rate  behavior,  a  more  realistic  description  of 
the  irrecoverable  deformation  should  also  account  for  coupling  between  plasticity, 
damage  and  viscous  effects.  A  discussion  of  other  choices  for  internal  variables  is 
given  in  this  section. 

We  denote  by  h  e  Kn  the  vector  of  n  scalar  internal  parameters  in  the  Eulerian 
description.  All  the  constitutive  functions  F,  a  and  0  will  depend  also  on  the  internal 
vector  h  and  the  constitutive  equation  (65)  becomes 


Wc  complete  the  model  (57)  by  providing  the  evolution  equation  for  the  internal 
variable  h .  Wc  assume  that  the  total  (particular)  derivative  of  h  is  a  function  of  h 
and  the  stretch  tensor,  i.e. 


(58) 


If  H  is  isotropic  with  respect  to  D  then  it.  depends  on  D  through  its  invariants*  For 
simplicity,  we  further  suppose  that  H  is  a  function  of  only  the  first  two  invariants 
of  D *  Le.  H  —  H(d,e.h). 
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A  first  choice  for  the  description  of  hardening  parameter  is  the  density  p .  In 
this  way,  we  capture  the  effect  of  pore  closing /col  lapse  on  the  response  of  geo¬ 
logical /cementitious  materials.  If  we  set  n  =  1  and  h  =  p  then  the  evolution 
equation  for  the  hardening  parameter  (66)  reduces  to  the  continuity  equation  (12), 
and  H  =  H  (d,  e,  p)  —  — pd .  In  many  applications,  the  volumetric  (plastic)  strain 
—  p/p{\  *-  1,  is  used  instead  of  the  density.  In  this  case  (66)  becomes 

(^  =  (^+v-Vn  =  -(l+fi)d,  (59) 

and  H  =  H(d,e,p)  =  —  (//  +  l)d. 

Other  authors  (c.g.  Cristcscu  and  Hunsche,  1998)  use  as  hardening  parameter 
the  visco-plastic  work  per  unit  volume  in,  i.e.  n  —  1 ,  h  =  uk  where  w  =  a  D 
From  (57)  wo  deduce 

—  =  —  +  v  ■  Vw  =  aid,  e,w)d  +  ft{d.  e,  w)e2,  (CiO) 

dt  at 

with  H  —  H(d*e,  w)  —  a(dre,w)d  4-  /3(d, e, w)e2. 

To  capture  the  combined  effect  of  plasticity  and  damage,  we  propose  the  following 
extension  of  the  revisited  Bingham  model  (68)- (29): 


r 

n1  - 

n  ,  *(M)" 

r  \d‘\  \ 

tj  — • 

,  k'i  < 

k{p,5) 

if  [D'|  ^0 
if  |Z>r|  =  0, 


(61) 


f  tr(cr)  f 3  =  — pc(/x)  +  (A  +  2r}/3)div{v)  if  d£u(i?)  <  0 
div(v)  <  0,  {  (62) 

[  tr(cr)/3  >  —  pc(p)  if  drii(v)  —  0, 

with  two  internal  parameters:  the  volumetric  strain  p  and  a  scalar  damage  parameter 
5,  The  damage  parameter  is  associated  to  the  loss  of  cohesive  strength  due  to  air 
pore  collapse.  Following  Holmquist  ct  al.  (1993),  we  consider  that  the  shear  flow 
stress  is  of  the  form 


=  a0\/2/3  mmlSjnax,  4(1  -  rf)  +  B(pc(p)/er0)N }.  (63) 


The  function  pc{fi),  is  plotted  in  Figure  1  (expression  given  in  section  6.1). 

The  evolution  equation  for  p  is  (59)  while  for  the  damage  parameter  we  consider 


dS 

dt 


68 

—  — -  +  v  ■  Vtf  = 

dt 


A  (fi)' 


(64) 


with  A(/t)  =  \/2/3  m  ax  {To,  D\[{p,.{p)  +  T)ja^\!>- } .  Note  that  these  evolution  equa¬ 
tions  can  be  written  in  the  form  (66)  with  n  —  2,  h  =  (ft,  5)  and  H(d,e.h)  = 
(-(1  +  n)d,c/A{fi)). 
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Figure  7:  Left:  the  evolution  of  the  internal  parameters  t  — *  (p(f),  <5(£))).  Right: 
the  evolution  of  the  stress  t  — ►  (— p(£)*9(0)  during  the  second  experiment  ( V/Lq  = 
100/sec)  for  the  model  with  damage  (dashed  line)  and  without  damage  {solid  line). 

We  further  illustrate  the  response  of  this  model  for  concrete  subject  to  uniaxial 
compression  conditions.  The  law  of  variation  of  pc(p)  is  taken  as  that  reported  by 
Holmquist  et  ah  (1993),  which  was  given  in  the  previous  section  (see  also  Figure  1). 
We  take  r0  ~  0,01,  Dy  —  0.04,  D2  —  l-  and  for  the  tension  cut-off  pressure  T  (see  51, 
we  take  the  value  reported  in  Holmquist  et  al.  (1993),  i.c,  T  —  4, MPa.  In  Figure  7 
(left)  is  shown  the  evolution  of  the  damage  as  a  function  of  the  compaction  level  ft 
for  the  loading  rate  V/Lq  =  100/s  as  in  the  second  experiment  but  with  ///,„  =  0.15. 
Note  that  during  loading,  damage  is  continuously  increasing  until  it  reaches  the  value 
Sjin  =  0.6,  Figure  7  (right)  shows  in  the  p  —  q  plane,  the  stress  variation  during 
the  same  test  {V/Lq  =  100/s)  in  the  ease  as  described  by  the  revisited  Bingham 
model  with  damage  (dashed  line)  and  without  damage  (solid  line).  As  expected, 
the  presence  of  damage  causes  a  decrease  in  the  yield  stress. 

2.8  Conclusions 

A  general  methodology  for  constructing  fluid-type  constitutive  equations  for  descrip¬ 
tion  of  the  combined  effects  of  plasticity,  damage  and  viscous  effects  on  the  behavior 
of  cementitious  and  particulate  media  when  subjected  to  large  deformations  and 
high  strain  rates  was  proposed.  The  clastic  deformations  were  neglected  and  the 
representation  of  the  dynamic  state  of  stress  in  the  material  was  obtained  from 
classical  yield  conditions  for  plastic  solids  using  two  methods:  (1)  the  viscoplas¬ 
tic  regularization  and  (2)  a  Maxwell  type  stress  superposition  method.  Examples  of 
constitutive  equations  obtained  from  classical  two-invariant  yield  conditions  F  using 
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both  procedures  arc  presented.  In  general,  even  if  we  stm  t  from  the  same  yield  con¬ 
dition,  the  models  obtained  using  these  two  methods  are  different.  In  the  case  of  the 
Mises  yield  condition  which  is  a  function  only  of  the  second  invariant,  using  cither 
model  we  obtain  the  classical  Bingham  fluid  model.  To  better  capture  the  behavior 
of  porous  geological  or  cementitious  materials  under  severe  loadings  (high-pressure, 
high  strain  rate  behavior),  extensions  of  the  Bingham  classical  model  were  proposed 
using  both  procedures.  Those  rigid  viscoplastic  fin  id  models  are  compressible,  their 
flow-no  flow  threshold  and  subsequent  flow  behavior  being  dependent  on  the  density 
or  compaction  level. 

We  showed  that  for  some  yield  conditions  F ,  certain  kinematic  restrictions  need 
to  be  satisfied.  In  the  case  of  the  Mises  yield  condition,  the  kinematic  constraint 
that  ought  to  be  fulfilled  is  the  incompressibility  condition.  However,  for  certain 
classical  yield  functions,  such  as  Drucker-Pragcr  the  resulting  kinematic  constraints 
are  physically  unacceptable  and  we  cannot  arrive  at  consistent  rigid-viscoplastic 
fluid  type  formulations.  For  other  yield  conditions,  such  as  the  Cap  model  with  a 
standard  Drucker-Pragcr  envelope,  the  derived  fluid-type  model  has  discontinuities 
with  respect  to  the  deformation  rate,  which  makes  the  model  inconsistent. 

We  illustrated  the  response  of  the  proposed  models  for  uniaxial  compression 
conditions  and  compared  their  respective  responses. 

Extensions  of  these  models  that  include  coupling  between  damage,  plastic  and 
viscous  effects  are  proposed  and  their  predictive  capabilities  and  limitations  were 
discussed. 
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3  Task  2:  Numerical  modeling  of  projectile  pene¬ 
tration  into  compressible  rigid  visco-plastic  me¬ 
dia 

3.1  Introduction 

One  of  the  many  challenges  associated  with  modeling  penetration  is  that  very  limited 
experimental  data  arc  available.  Although  progress  has  been  made  in  instrumenta¬ 
tion  and  diagnostic  techniques  that  allow  accurate  measurements  of  impact  velocity 
and  projectile  deceleration,  in-target  instrumentation  is  limited.  Thus,  information 
concerning  the  stress-time  response  in  the  target  during  penetration  is  lacking  (see 
for  example,  Gran  and  Frew  [lj).  To  gain  insights  and  fundamental  understanding 
of  the  phenomena  involved,  numerical  simulations  of  penetration  processes  play  a 
key  role. 

For  modeling  the  deformation  of  the  target  during  penetration  Lagrangian  and 
Eulcrian  methods  have  been  used  {see  for  e.g.  the  survey  by  Hertel  [2]).  Each  of 
these  methods  has  advantages  and  disadvantages  as  a  tool  to  model  penetration. 
With  Lagrangian  methods,  the  conservation  equations  are  solved  using  either  finite 
difference  or  finite  element  techniques  on  a  mesh  which  moves  with  the  material. 
The  material  interfaces  arc  tracked  in  a  sharp  fashion.  However,  severe  deformation 
leads  to  mesh  tangling  or  distortion,  hence  considerable  complexity  is  enjoined  by 
the  need  for  remeshing.  Mesh- loss  methods  (e.g.  [3]-[6] )  or  a  combination  of  finite 
element  methods  with  embedded  boundary  tracking  and  local  enrichment  (c.g.[7j  ) 
have  emerged  in  recent  years.  When  using  Eulcrian  methods,  the  conservation  equa¬ 
tions  are  solved  using  either  finite  difference  or  finite  volume  techniques.  However, 
when  using  such  methods  the  main  challenge  is  to  be  able  to  track  moving  material 
boundaries  (i.e.  problems  related  to  recognition  of  material  interfaces  or  surfaces). 
Such  difficulties  arc  the  driving  force  for  the  development  of  Arbitrary  Lagrangian 
Eulcrian  (ALE)  and  level  set  methods  (see  for  e.g.  [2,  8]). 

It  is  to  be  noted  that  the  material  models  used  in  conjunction  with  the  aforemen¬ 
tioned  methods  apply  mainly  to  penetration  into  metallic  materials.  It  was  shown 
(sec  for  instance  (9,  10])  that  the  use  of  an  Eulcrian  setting  and  of  a  fluid- type 
constitutive  equations  presents  clear  advantages  over  other  methods  employed  for 
modeling  penetration  in  metallic  materials.  However,  implicit  in  these  fluid-type 
models  is  the  hypothesis  of  incompressibility,  thus  these  models  cannot  describe 
the  irreversible  volumetric  changes  observed  in  geologic  or  cementitious  materials. 
Indeed,  from  post- test  observations  of  the  density  distribution  in  the  target  it  can 
be  clearly  seen  that  when  subjected  to  dynamic  loading,  porous  and  particulate 
materials  show  large  irreversible  volumetric  and  shear  deformations,  the  behavior 
being  highly  sensitive  to  the  strain  rate  (e.g.  [11,  12,  18]).  Batra  and  Gobinath  [14] 
have  proposed  a  compressible  model  that  accounts  for  the  rate  dependence  of  the 
shear  response.  Hardening  effects  arc  neglected  and  the  pressure-density  response  is 


33 


given  by  a  gas  type  state  equation  which  cannot  capture  plastic  volumetric  effects. 
Recently ,  for  modeling  penetration  in  concrete,  a  purely  hydrodynamic  treatment 
was  used  in  [15]  while  LS-DYNA  finite-element  simulations  have  been  reported  by 
[16].  Under  this  grant,  support,  computational  techniques  for  modeling  steady-state 
flow  of  a  compressible  rigid  viscoplastic  fluid  have  been  developed.  Illustration  of 
the  application  of  these  methods  was  provided  for  cementitious  materials  considered 
to  obey  a  revisited  Bingham  type  model  [18],  which  accounts  for  plastic  volumetric 
deformations.  The  reader  is  refered  to  the  journal  paper  Cazacu  ct  al  [17]  in  which 
this  work  is  presented  in  detail. 

In  the  following  we  will  present  the  computational  methods  developed  for  model¬ 
ing  non-steady  penetration  of  a  rigid  projectile  into  cementitious  targets.  To  capture 
the  solid-fluid  transition  iri  behavior  at  high  strain  rates  observed  in  such  media  and 
account  for  damagc/plastieity  couplings  and  how  these  dissipative  mechanisms  are 
influenced  by  the  strain  rate,  we  consider  a  compressible  rigid  viscoplastic  fluid 
model,  that  we  called  Bingham  revisited  {see  section  2,4.2  of  this  report  and  the 
journal  paper  in  which  we  published  the  model  [18]), 

A  hybrid  time-discretization  is  used  to  model  the  non-stationary  flow  of  the  the 
target  material  and  the  projectile- target  interaction,  i.e.  an  explicit  Euler  method 
for  the  projectile  equation  and  a  forward  (implicit)  method  for  the  target  boundary 
value  problem.  At  each  time  step,  a  mixed  finite-element  and  finite-volume  strategy 
is  used  to  solve  the  boundary  value  problem  in  the  target.  Specifically,  algorithms 
for  solving  the  nonlinear  variational  problem  involving  the  velocity  field  via  the  finite 
clement  approximation  are  developed  while  finite-volume  techniques  are  adopted  for 
solving  the  hyperbolic  mass  conservation  and  damage  evolution  equation. 

ft  is  to  be  noted  that  when  using  a  rigid  compressible  viscoplastic  fluid  equa¬ 
tion,  additional  difficulties  arise  due  to  the  non-smoothness  of  the  plastic  terms. 
Moreover,  since  in  the  constitutive  description  for  the  target  material  a  rigid  un¬ 
load  ing  hypothesis  is  considered,  the  variational  inequality  to  be  solved  has  unilateral 
constraints  for  the  velocity  divergence.  Thus,  we  cannot  simply  make  use  of  finite 
element  techniques  developed  for  Navicr-S takes  fluids.  Therefore,  to  solve  the  veloc¬ 
ity  problem,  a  decomposition-coordination  formulation  coupled  with  the  augmented 
lagrangian  method,  developed  by  [19,  20]  (for  a  comprehensive  review  see  also  [21]  ) 
for  the  Bingham  fluid  is  adapted  here.  The  velocity  problem  and  the  methods  used 
to  solve  it  are  similar  to  those  developed  in  |17]  for  modeling  steady-state  flow  of  a 
revisited  Bingham  fluid. 

It  is  shown  that  using  the  material  model  and  algorithms  presented  in  this  paper, 
it  is  possible  to  provide  details  concerning  the  stress  and  kinematic  fields  in  the 
target.  Specifically,  the  density  changes  around  the  penetration  tunnel,  the  shape 
and  location  of  the  rigid /pi  as  tic  boundary,  the  extent  of  damage  due  to  air  void 
collapse  in  such  materials  arc  accurately  predicted.  Moreover,  the  numerical  model 
is  robust  and  computationally  inexpensive.  This  is  due  to  the  fact  that  in  modeling 
penetration  using  a  rigid- type  model  only  a  small  zone  around  the  projectile  is 
deformed  and  thus  we  can  restrict  the  computations  to  a  small  bounded  domain 
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around  the  projectile.  Another  advantage  of  using  a  fluid- type  approach  is  that  the 
computations  arc  done  in  a  fixed  domain  and  there  is  no  need  to  mesh  different 
domains  at  each  time  interval. 

One  of  the  many  challenges  associated  with  simulating  penetration  in  cementi¬ 
tious  materials  is  related  to  the  lack  of  experimental  data  for  the  high  pressures  and 
high  strain  rates  involved  (sec  e.g.[22].  To  assess  the  effects  of  the  constants  involved 
in  the  material  model  that  cannot  be  determined  explicitly  from  experimental  data, 
a  parametric  study  was  also  performed.  Thus,  we  analyzed  the  sensitivity  of  the 
penetration  depth  calculations  with  respect  to  shear  and  volumetric  viscosities,  cut¬ 
off  yield  limit,  and  frictional  constants.  The  results  of  this  study  provides  insights 
into  the  relative  importance  of  the  plastic  and  viscous  effects  on  the  penetration 
process. 


3.2  Constitutive  model  for  the  target  material 

In  this  section,  we  give  an  overview  of  the  material  model  used  for  the  target  . 
To  capture  the  solid-fluid  transition  in  behavior  at  high  strain  rates  observed  in 
cementitious  porous  media,  account  for  damagc/plasticity  couplings  and  how  these 
dissipative  mechanisms  are  influenced  by  the  strain  rate  we  consider  a  compressible 
rigid  viscoplastic  fluid  constitutive  equation.  The  general  form,  proposed  in  [18] (see 
also  section  4  of  the  present  document)  is: 


D(v)  e  C 


■{ 


<T  =  q(D(v),  h)I  +  /?(£?(«).  h)D'(v) 
F(rr,h)  <  t) 


if  |D(tOM  0, 
if|D(v)|  =  0, 


(65) 


In  (65),  a  denotes  the  Cauchy  stress  tensor,  D  is  the  rate  of  deformation  tensor, 
D  =  D(v )  —  (Vu+VJ  v)/‘2  ( v  denotes  the  velocity  field)  w-hile  D'  —  D  —  t.r(D)/3I 
is  its  de viator.  I  is  the  identity  tensor,  h  g  R"  is  the  vector  of  n  scalar  internal 
variables.  F[<j.h)  is  the  yield  function  and  C  C  R|*3  represents  a  set  of  kinematic 
constraints,  while  a  and  ,3  are  scalar  functions  of  their  arguments.  The  invariance 
under  superposed  rigid  body  motions  requires  that  u  and  0  are  isotropic  functions 
of  D 

The  presentation  of  the  model  (65)  is  completed  by  providing  the  evolution 
equations  for  the  internal  variables  h.  It  is  assumed  that  the  total  (particular) 
derivative  of  h.  is  a  function  of  h  and  the  stretch  tensor  D ,  i.e. 


vh-at  dm.  m 

If  there  is  flow  (i.c.  in  the  visco- plastic  regime  where  D  0),  then  F(a(JD, 
fj(D,  h)D\h )  >  0,  for  all  D  G  C\0*  It  should  be  noted  that  at  difference  with  a 
classic  fluid  constitutive  equation,  for  a  rigid  viscoplastic  material  the  constitutive 
functions  a  and  fl  are  not  defined  and  cannot  be  prolonged  by  continuity  in  D  —  0. 
To  ensure  continuous  transition  between  flow  and  no- flow  states,  it  is  required  that 
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<•»(■.  h),  ft{-.  h)  axe  continuous  on  C\0.  The  continuity  of  the  transition  between  flow 
and  no  flow  reads:  F(a(D,h)I +  ft(D,h)D'.h) —*  t),  for  D  — »  0  with  D  €  C\0, 
while  the  fact  that  the  dissipated  power  has  to  be  nonnegative  leads  to  a  D  = 
a(D ,  h)tr(D)  +  ft(D.h)\D'\2  >  0  for  all  D  £  C. 

Two  general  procedures  for  determining  the  expressions  of  the  constitutive  func¬ 
tions  F,  a,  and  ft  such  that  the  model  (65)- (66)  is  consistent  (i.e.  the  above  com¬ 
patibility  conditions  arc  satisfied)  are  provided  in  [18].  Although  all  the  numerical 
methods  and  techniques  developed  in  this  paper  arc  valid  for  the  constitutive  equa¬ 
tion  (G5)-(66)  in  its  general  form,  in  this  paper  we  consider  a  specific  compressible 
rigid  viscoplastic  model.  In  this  model  it  is  assumed  that  yielding  under  deviatoric 
conditions  is  of  Von  Miscs  type  will)  yield  limit  k  depending  on  two  internal  vari¬ 
ables:  the  current  density  p  and  the  damage  parameter  (5,  i.e.  n  =  2,  h  =  (p;  <5) . 
Thus,  the  effect  of  the  density  on  flow  is  taken  into  account  while  the  damage  pa¬ 
rameter  is  associated  to  the  loss  of  cohesive  strength  due  to  air  pore  collapse. 

For  hydrostatic  conditions,  the  response  is  described  by  a  pressure-density  re¬ 
lationship  of  the  form:  p  —  pr(p ).  Note  that  a  full  equation  of  state  which  gives 
a  complete  thermo-dynamical  characterization  is  not  available  for  porous  cemen¬ 
titious  materials  (see  for  example  the  very  recent  survey  and  experimental  study 
on  concrete  [22]).  This  is  due  to  the  fact  that  using  the  current  testing  methods 
(detonation  and  flyer-impact  testing)  only  data  pairs  of  pressure-density  points  can 
be  obtained. 

Data  available  for  cementitious  materials  (see  for  instance  [23])indicate  a  stiff 
response  upon  unloading.  Thus,  the  reversible  decrease  of  volume  can  be  neglected 
anil  the  unloading  process  can  be  considered  rigid  (the  density  acquired  at  the  end 
of  the  loading  process  is  preserved  during  unloading).  Hence,  the  yield  condition  is 
of  the  form: 

F{p, i i, p, S)  =  ^\^q-  k(p,«5)15.  +  [p  +  pr(p)]“.  (67) 

where  [s]+  =  (|ar|  +  x)/2,  [x]_  =  (|x|  -  x)/2  and  p  -  ~  ^  ■,  r/  -  |er'|.  Using 

this  yield  function  and  the  superposition  technique  (he.  compressible  Navicr-Stokes 
fluid  superposed  on  an  associated  plastic  solid  described  by  (67))  we  can  deduce,  (see 
[18]),  the  following  constitutive  equation  called  the  revisited  Bingham  fluid  model. 
The  deviatoric  part  of  this  model  is: 


{ 


k'l< 


2  v(p)  + 

K(pJ) 


«(M)' 
\D'\  . 


D1 


if  \D‘\  ?() 
if  \D'\  =  0, 


(68) 


while  its  spherical  part,  is: 

/r(er)/3  =  -Pr(p)  +  IHp)  +  2//{/?)/3]  div(v)  if  div(u)  <  0 


div(u)  <  0, 


tr(er)/ 3  >  -p,:(p) 


if  div(u)  —  0. 


(69) 
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Figure  8:  Representation  of  the  domain  T>o  for  the  tunnelling  stage  of  penetration. 


Note  that  the  classical  Bingham  fluid  is  recovered  if  in  (68)  the  incompressibility 
condition  (he.  divv  =  0)  is  imposed. 

The  kinematic  restriction  is  C  =  {I?  ;  tr(D)  <  0}  i,e.  the  flow  is  compress¬ 
ible  and  unloading  is  rigid.  In  many  solid-type  models,  pc  and  k  are  expressed  as 
functions  of  the  compaction  level  (or  volumetric  strain)  p  =  p/pQ  -  1  instead  of  the 
density  p  (here  po  is  the  initial  density  prior  to  loading).  Note  that  (68)-(69)  can  be 
recast  in  the  form  that  is  generally  used  in  hydrocode  calculations  by  eliminating  // 
between  k  =  h:(/t,  <S)  and  pr  —  pc(/x). 


3.3  Geometry  and  fields  equations 

When  a  penetrator  impacts  a  semi- infinite  geological,  or  cementitious  target,  the 
material  is  displaced  radially  and  as  a  result  a  tunnel- shaped  crater  is  formed  in  the 
target.  The  entrance  portion  of  the  crater  is  generally  produced  by  spallation  in  the 
vicinity  of  the  impact  site.  The  crater  entrance  is  typically  wide  and  shallow,  but  it 
rapidly  evolves  into  a  circular  tunnel. 

Our  objective  is  to  model  the  tunnelling  phase  of  the  penetration,  when  the  full 
length  of  the  projectile  is  in  the  target  and  the  target  material  flows  around  the 
projectile.  It  is  assumed  that  the  projectile  is  rigid. 

Let  us  denote  by  Vq  and  %  the  domains  occupied  by  the  projectile  and  by  the 
semi-infinite  tunnel  in  the  coordinates  Oxyz  which  arc  linked  to  the  projectile  (see 
Figure  i  ). 

Let  be  the  coordinates  attached  to  the  target  and  let  —V{t)e$  be  the 

velocity  of  the  projectile  at  the  moment  t  >  0.  The  two  system  of  coordinates 
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coincide  at  /  =  0  and  we  have 


X\  =  X,  x2 


—  V,  x:i  =  2  -  f 
J  0 


V($)  ds. 


The  target  occupies  the  variable  domain  Dt  which  is  the  whole  space  R:i  without  the 
projectile  Vt  —  {fai}x2lx^);  (xi,X2,x3  +  /  V'(s)  ds)  £  A}  and  the  infinite  tunnel 

J{ ) 


%  —  {  fa  1 1  x2 ix$)}  fa  1 t  x2 1  ^'3  + 


f  V(s)  ds)  £  %}  behind  it,  i,e.  A  =  R*  \  (Vt  U7J), 
Jo 

and  it  corresponds  to  2>0  =  R*\(Au7[i)  in  the  coordinates  linked  to  the  projectile* 
We  denote  by  T]  F^  the  part  of  the  boundaries  of  Tt  and  Vt  which  are  in  contact 
with  the  target,  lc.  &Dt  =  r\  U  Ff>,  which  corresponds  to  dt>0  ~  A  U  F2,  in  Oxyz 
coordinates. 

Next,  wc  present  the  equations  governing  the  non  stationary  flow  of  a  rigid 
viscoplastic  material  described  by  the  constitutive  equations  (67)and  £68)- (69).  In 
order  to  write  the  conservation  laws  in  the  fixed  domain  At  instead  of  the  variable 
domain  A,  we  denote  by  v  :  [0,T|  x  Aj  R;i  the  relative  velocity  field  given  by 

v(t,z,yyz)  =  v(t,xux2,x^)  +  V(t)ex, 

where  v  is  the  Eulcriao  velocity  field  in  A<  The  other  fields,  i.e,  the  stress  <r.  density 
p  mid  damage  5  will  have  the  same  notations  in  both  systems  of  coordinates. 

The  momentum  balance  law  in  the  Eulcrian  coordinates  Oxyz  reads 

p(/)[0fi;(O  +  (v(£)  *  V)v(t)  —  K(t)eJ  -  div*r(£)  =  0  in  Ah  (70) 

the  continuity  equation  is 

"  +  v  ■  Vp  =  -pdiv{v),  in  T>0m  (71) 

at 

while  the  damage  evolution  law  can  be  written  as 

-+v-V6  =  H(p,6.D(v))  in  T>0.  (72) 

at 

The  conditions  at  infinity  for  the’  target  in  the  coordinates  linked  to  the  projectile 
read: 

v{tyx,y,z)  =  V{t)em  for  ||(ar,y,^)||  -» +oo,  (73) 

p(t,x,y,z)  =  p°,  S(t,x>y,z)=  0  for  z  —  -oo.  (74) 

The  mechanical  model  considered  for  the  target  material  is  rigid- viscoplastic. 
Thus,  with  the  exception  of  a  domain  around  the  projectile  where  viscoplastic  de¬ 
formations  occur,  the  target  material  is  in  rigid  motion.  Hence,  the  domain  Ah 
affected  by  the  impact  event  ,  can  be  restricted  from  the  whole  space  to  a  bounded 
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region  (depicted  in  Figure  1  with  interrupted  lines).  To  enable  reasonable  computa¬ 
tional  effort  and  yet  ensure  that  we  capture  the  viscoplastic  region  in  its  entirety  (i.c. 
the  boundary  conditions  at  infinity  arc  accurately  described),  we  limit  the  extent 
of  this  domain  to  5  projectile  radii.  This  assumption  appears  to  be  supported  by 
experimental  evidence.  Indeed,  post-test  observations  indicate  that  the  tunnel  is  of 
the  order  of  the  projectile  diameter.  For  both  grout  and  concrete  targets  a  change 
in  density  in  the  region  around  the  penetration  tunnel  radially  outward  from  the 
edge  of  the  tunnel  to  a  distance  of  1  -  1.5  projectile  diameters  was  reported  (see 
Jones  ct  al  [24] ) . 

Using  this  assumption  the  boundary  conditions  at  infinity  (73)-(74)  are  replaced 
by 

«(i)  =  V(f)e*  on  r0l  p(t)  =  p°,6(t)  =  Oon  rin.  (75) 

On  the  surface  in  contact  with  the  tumid,  F  i,  the  stress  vector  vanishes 


a{t ,  x)n  =  0, 


(76) 


where  n  stands  for  the  outward  unit  normal  on  dV0,  while  frictional  contact  its 
assumed  on  the  boundary  of  the  projectile  F^- 

The  mechanics  of  friction  at  high  speed  of  sliding  is  very  complex.  Much  of  the 
work  reported  in  the  literature  is  at  lower  speeds  or  pressures  than  those  occurring 
during  a  penetration  event.  In  view  of  this,  in  our  analysis  we  assume  that  a 
Coulomb  friction  law  applies.  However,  in  order  to  capture  the  rate  effect  on  the 
frictional  contact  between  the  concrete  target  and  the  metallic  penetrator,  a  slip 
rate  dependent  friction  coefficient  is  assumed,  he.: 


Wn  =  0, 


'  v,  ~  0  =>  Iff, |  <  /loj, 

V,  ^  0  =>  CT,  =  ~T\an\~r 

[Vf| 


(77) 


where  uu  —  v  ■  n  is  normal  velocity,  an  —  an  *  n  is  the  normal  stress,  u*  =  v  —  unn 
is  the  tangential  velocity,  er*  =  a  —  [a  *  n)n  stands  for  the  tangential  stress  and 
JF  is  the  friction  coefficient.  According  to  (77),  the  tangential  (friction)  stress  is 
bounded  by  the  normal  stress  multiplied  by  the  friction  coefficient  T .  If  such  a 
limit  is  not  attained  sliding  cannot  occur;  otherwise  the  friction  stress  is  opposite  to 
the  slip  rate.  The  friction  coefficient  will  also  be  considered  variable.  Experimental 
observations  indicate  that  the  friction  coefficient  depends  on  the  slip  rate  \vt\  i.c. 
T  =  JF(|u,|).  The  simplest  law  of  variation  of  JF  on  the  slip  rate  is  a  discontinuous 
jump  from  a  “static”  value  (for  |u,|  =  0}  down  to  a  “dynamic"  or  “kinetic”  value 
(for  |rjf]  ^  0).  In  this  work,  we  will  consider  a  smooth  and  decreasing  function 
T  —  I)  of  the  slip  rate. 

If  wc  denote  by  M  the  mass  of  the  projectile  then  the  fluid  structure  interaction 
equation  is 

MV(t)  =  -  f  cr(t)n  *  e*.  (78) 

Jr2 
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Filially,  the  initial  conditions  arc 


v(0)  =  vn,  p(0)=Po,  «(0)  =  du,  V(0)=V0  (79) 

where  Vq,Pq,  and  £()  are  the  initial  relative  velocity,  density,  and  damage,  respec¬ 
tively.  We  suppose  that  the  initial  fields  are  smooth  and  compatible  with  the  bound¬ 
ary  conditions,  i.e.  no  discontinuities  are  generated  initially. 

3.4  The  algorithm 

3.4.1  Time  discretization 

Let  At  he  the  time  step  and  let  us  denote  by  vk.  erk,  pk,  Sk,  Vk  the  values  of 
v{kAt),  tr(kAt),  p(kAt),  S(kAt)  and  V(kAt).  Suppose  that  we  have  computed  all 
these  values  at  t  =  kAt.  The  time  explicit  (forward)  Euler  scheme  for  the  projectile 
equation  (78)  will  then  provide  Vk+1 

Vk+i  =  yk  _  &  f  akn  .  (80) 

M  Jrl 

and  the  time  implicit  (bachward)  Euler  scheme  for  (70)-  (72)  gives  the  following 
nonlinear  equations  for  el'+l .  crA+1 ,  pk+l  and  5k+l 

/jk+1[«*+i  +  A tvk+'  ■  Vuk+I)  -  At  div  <rk+l  =  pk+'fk  in  Va,  (81) 

pk+l  +  At  div(/3k+1uk+1 )  —  pk  in  T>0,  (82) 

Sk+l  +  AtV6k+1  ■  ufc+1  =  6k  +  A t.H(pk+l,6k+l,D(vk+l))  in  V0.  (83) 

where  fk  —  vk  +  ( Vk+]  —  Vk)et.  To  this  set  of  equations,  we  have  to  add  the 
constitutive  equations  (G7),(G8)-(69)  which  link  all  three  unknowns  i.e.  t>k+l,  rr*'  +  11 
pk+l.  The  conditions  at  the  "artificial”  boundary  of  P(J.  which  replace  the  conditions 
at  infinity,  arc 

vk+ 1  -  V*+lee  on  T0,  pk+l  =  p°,  Sk+l  =  0,  on  Tm.  (84) 

Moreover,  uk+l.  u*+I,  a*+1  and  crk+i  have  to  satisfy  the  friction  law  (77). 

Setting 

Wfc+1  =  {v  €  H'CDq)3  ;  v  =  K*+leE  on  T0,  =0  on  T2}, 

the  kinematic  admissible  set  is: 

Vk+1  =  {ti£  Wk+l  ;  divv  <  0  in  I>d}, 
and  the  variational  formulation  for  the  velocity  field  i/'+1  G  V*+1  is: 

f  pk+lvk+l  ■  (v  -  u*+1)  +  A t{[  pk+'vk+1  ■  Vvk+l  ■  (u  -  t>*+‘)+ 

Jv o  Jt> o 
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/  2 r)(pk+l)D'(vk+l):(D,(v)~D'(vk+l))+  f  K(pk+i, d*+1)(|Z)'(«)j  -  |DV+1)I) 

Jt>0  ^Pr> 

+  f  [(A(p*+1)  +  ^J7(p*+1))div«fc+1  -  pc(pfc+l)]div(tJ  -  vAf+l) 

Jv „  3 

+  /  ^(|«i+,|)k:+,l<W  -  |„J«|)}  >  /  p‘+‘A  ■  (v  -  v"*'),  (85) 

j  Vi  JVa 

for  all  v  €  V*+1. 

3.4.2  The  general  scheme 

Lot  us  lix  the  iteration  in  time  k.  The  numerical  algorithm  for  solving  the  nonlinear 
equations  (8 1 )-{ 83) ,  consists  of  solving  alternatively  the  variational  inequality  (85) 
for  the  velocity  field  and  the  equations  (82)-(83)  for  the  density  and  damage  fields. 
More  precisely,  we  will  distinguish  two  problems  :  the  ” velocity  problem”  and  the 
’'density/damage  problem”. 

For  the  velocity  problem  (at  iteration  n  of  the  algorithm),  we  assume  that  pk+l  = 
pk+i,n  tmfj  £*+i  _  £*+i,n  jy-g  gjVCn!  i.c.  the  distributions  of  p  =  pf(pk+l),V  — 
7l(pk+ 1  )i  A  —  A(pfc+I)  and  k  =  n(pk+i ,  6t+1)  are  known,  and  we  find  the  velocity 
o;'+I  —  vfc+l,n+1  €  V  by  solving  the  variational  inequality  (85). 

The  density /damage  problem  (at  iteration  n  of  the  algorithm)  consists  in  finding 
the  density  and  damage  fields  i.c.  pfc+1  =  pfc+1'n+1  and  <Sfc+1  =  <5fc+l'n+1  solutions  of 
(82)  and 

6k+l  +  Atdiv(6k+lvk+1)  =  Hkn  in  VQ.  (86) 

where  Hk  n  =  6k  +  Al\H(p/c+l'n ,Sk+],n ,  D(vk+l))  +  <5fc+1,n  div  vfc+1],  assuming  that 
vfc+1  =  jg  known. 

Since  we  will  suppose  that  the  time  iteration  is  fixed,  in  the  following  we  will 
omit  the  subscript  k  +  1.  Moreover,  in  order  to  simplify  the  notations,  we  will  also 
omit  the  subscript  n  +  1. 

3.4.3  The  velocity  problem 

It  is  to  be  noted  that  in  the  ease  of  a  rigid  compressible  visco plastic  fluid,  additional 
difficulties  arise  from  the  non  differentiability  of  the  plastic  terms  and  from  the 
unilateral  constraints  on  the  velocity  divergence.  That  means  that  we  cannot  simply 
make  use  of  the  finite  element  techniques  developed  for  Navicr-Stokcs  fluids  (see  for 
instance  [25,  26]).  To  solve  the  velocity  problem,  a  decomposition-coordination 
formulation  coupled  with  the  augmented  lagrangian  method  (sec  [19,  20])  will  be 
adapted  for  the  material  model  considered  (67), (68)- (69).  After  time  discretization, 
the  variational  formulation  in  terms  of  velocities  is  similar  to  that  corresponding 
to  the  stationary  flow  of  a  revisited  Bingham  fluid  given  in  [17].  Although  the 
numerical  techniques  developed  here  for  the  velocity  problem  are  very  close  to  that 
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proposed  in  [17],  in  the  following  we  present  its  main  features  and  refer  the  reader 
to  [17]  for  further  details. 

Do  is  discretized  by  using  a  family  of  triangulations  (Th  )h  marie  of  finite  elements 
(of  degree  2),  where  /;  >  0  is  the  discretization  parameter  representing  the  greatest, 
diameter  of  a  triangle  in  Th.  The  finite  element  space  W>,,  which  is  an  internal 
approximation  of  W  will  be  denoted  by  W4,  while  V/,  =  Vfl  WV  Then,  the  velocity 
problem  is  discretized  by  considering  V/t  €  Vj,  which  satisfies  the  variat  ional  inequal¬ 
ity  (85)  for  all  v  €  V;,.  In  order  to  simplify  the  notations,  we  shall  omit  through) 
this  subsection  the  index  h. 

For  all  w  6  V/,  wo  define  Jw  :  V),  -*  R  as 

Jw{v)=  f  £|u|2  +  Af{/  J7|D'(v)|2+  /  +  ^Vdivu)2  +  f  «|D'(tO|  + 

J  -D0  *  J  Vo  JV0  \z  *V  JVq 

f  F(\W'\)\(T„(w)\\v,\+  I  p(w -V)w -v  -  f  pdivt?}  —  f  ()f  V.  (87) 
J  r3  Jv  Jv0  Jvq 

and  rewrite  the  velocity  problem  as  a  quasi- variational  problem 

v  €  14,  J„(v)  =  min  Jv(v),  (88) 

To  reduce  the  quasi- variational  problem  to  a  sequence  of  minimization  problems, 
we  use  the  following  iterative  algorithm  {I A)  :  given  vJ~{  €  Vtt  find  vJ  €  Vft,  the 
solution  of  the  following  minimization  problem  for  Jj{v)  = 

vJ  €  Jj(vJ)  —  min  Jj{v).  (89) 

veVf, 

Since  Jj  is  not  Gateaux  differentiable,  one  can  use  a  regularization  tcchnique(scc 
for  instance  [27,  28|).  With  this  technique,  the  material  is  not  completely  rigid 
anymore,  and  so  it  is  difficult  to  capture  accurately  the  shape  of  the  rigid  zone. 
To  overcome  this  difficulty  and  describe  the  rigid/ viscoplastie  boundary  sharply, 
we  adapt  the  decomposition-coordination  formulation  coupled  with  the  augmented 
lagrangian  method  originally  introduced  for  the  incompressible  Bingham  fluid  in  [19, 
20].  On  the  other  hand,  for  the  frictional  terms,  which  are  also  non-differentiable,  we 
shall  use  a  classical  regularization  technique  because  in  the  applications  considered 
(penetration  into  geological  materials)  the  slip  rate  is  always  non- zero. 

Let  us  consider  the  spaces  of  lagrangian  multipliers:  A*  ,  which  denote  tensor 
valued  polynomial  functions  of  degree  one  on  each  triangle,  and  0*  which  denote  real 
valued  polynomial  functions  of  degree  one,  respectively  .  The  augmented  lagrangian 
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Cj  :  Vh  x  A/,  x  0h  x  A/,  — *  R  is  : 

Cj(v.(T<Ot-i)  =  f  ~\v\2  +  At{  f  vfo\'2  +  f  (^  +  ^V2+  /  «M~ 

Jv0  z  Jt>q  ^  JVi I  \*  J/  Jpn 

[  pO  +  f  Jr(lW~1|)K(«J~1)|—7^  = — +  [  p{v3~l  •  V)uJ_1  -V  + 

Jv0  Jr,  VW-IP  +  «*  Jv° 

I  {D'(v)  —  7)  :  <r'  +  j  (0  -  div  v)  tracc(<r)}  -  /  pf  ■  v  + 

JVq  J  Do  0 

ro/jJD'H-Tl2  +  rH  |  6-  divu|2. 

(90) 

whore  r},  and  rD  are  strictly  positive  constants.  Here  7  e  Ah  stands  for  D'(v).  0  € 
0;,  stands  for  div  v.  and  cr  £  A;,  is  a  iagrangian  multiplier.  Thus,  the  minimization 
problem  (89)  becomes:  find  vJ  £  V*,,  erJ  €  Ah,  01  £  0?,  and  jJ  €  A,,  such  that 

sup  Cj(vJ,tr,8,y)  <  Cj(vk,<ri,(P,‘yi)  <  inf  £J(v,(ri1^,7J)>  (91) 

0,t r.T  *«V* 


In  order  to  solve  the  above  saddle  point  problem,  we  use  an  Uzawa-type  algorithm 
(UA).  The  interest  in  using  this  algorithm  is  that  it  transforms  the  non-diffcrcntiablc 
problem  into  a  sequence  of  completely  standard  computations. 

Let  first  define  the  functions  ///  and  /d,  useful  in  the  description  of  the  algorithm: 


fn(cr,a)  = 


1 


1  - 


«(p. 


2 T}(p)  +  a  l  \<t’\  . 


tr\  /w(cr ,  a)  =  - 


3p(p)  +  tr  a 


,3A(p)  +  2i]  (p)  +  a 


,  (92) 


where  [.r]+  =  (x  +  |x|)/2  and  [x]_  =  (|i|  — 1)/2  are  the  positive  and  the  negative  part 
of  [j],  respectively.  The  (UA)  algorithm,  starts  by  choosing  0^  1  =  tr-'-1,  1  —  0J~i 

and  7^_l  =  7-'-1  as  initial  guesses.  At  the  iteration  j,  rrU’.Tf.I,1, 0\Z\  are  already 
computed  and  we  look  for  vf  ~ 1 ,  <rJ~ 1 , 7^ " 1 ,  ‘ 1 . 

First,  we  compute  v,~ 1  solving  the  following  quadratic  minimization  problem 


v 


e  Wh, 


£*(* j  I . Tf-I )  =  tof 


Then,  we  compute  explicitly  7^  1  and  1  through  the  functions  ///  and  /„ 
7>_l  -  +  2rDl?(u|-1),2r0)>  Sf*  =  /„(<,*  +  2r„i>(tr|- ')), 2r„) 

and  we  update  er(_l  through  the  following  formulas: 

«*)'  =  (o’Ci1)'  +  2td{D‘{v3~')  —  7j-1), 

traced"1)  —  trace (o^Ij)  +  2rw(divuf~‘  -  0f“l). 

The  (UA)  algorithm  stops  at  i  =  when  the  error  e  =  II/.  iph)  <g 

ll1^  lka(fD) 

small  enough.  Then  we  update:  ~  1  —  cr(,“ 1 , 0J  =  0f,_I  and  7J  —  7^  l.  In 
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the  numerical  computations*  when  the  solution  v^~2  is  dose  to  t^-1  (for  instance 
when  At  is  small  enough),  a  single  pass  of  Uzawa-type  algorithm  (i.c,  =  1}  can 

be  enough  to  satisfy  the  stop  criterion  in  i.  This  fact  reduces  significantly  the  total 
computation  time. 


The  algorithm  (1A)  stops  at  j  =  jcT  when  the  error 


l|t)J  -tJ^'II^CP,,) 


is  small 


enough. 


3*4*4  The  density /damage  problem 

We  will  use  a  finite  volume  method  (sec  for  instance  (29])  to  discretize  the  continuity 
equation  (82)  and  the  damage  evolution  equation  (86)*  Let  denote  by  /C/,  the  finite 
volume  mesh,  which  is  given  by  a  family  of  disjoint  polygonal  connected  subsets  of 
!R*  such  that  V  is  the  union  of  the  closure  of  the  elements  of  fCh\  h  is  the  greatest 
diameter  o[  a  control  volume  in  fCh,  The  finite  volume  mesh  Kh  may  coincide  (or 
not)  with  the  finite  element  triangulation  Th.  In  the  simulations  presented  in  the 
next  section  fCfl  =  Th . 

The  finite  volume  discrete  space  is  the  space  of  piecewise  constant  functions,  i.e. 
we  arc  looking  for  the  solution  ph  of  (82)  and  4  of  (86)  as  {phK,  A  €  AC/,}  and 
{&hk\  A  €  /C/,},  respectively.  In  order  to  simplify  the  notations,  throughout,  this 
subsection  we  omit  all  the  indexes  h. 

Lei  us  consider  two  control  volumes  h\  P  6  A4  with  a  common  interface  Ikl  — 
Kf]P.  Let  n  be  the  unit  normal  vector  to  Ik p  oriented  from  A  to  P.  Then,  we 
define  the  flux  F(h\  P)  at  the  interface  I kp  as 

F(K,P)=  I  [v-n]+. 

'I  Ikp 

Note  that  at  least  one  of  the  two  fluxes  F(K,P)  and  F ( P.  A" )  is  vanishing.  If  we 
denote  by  M{h  )  the  set  of  all  neighbors  of  the  control  volume  A’,  then  the  finite 
volume  numerical  scheme  for  (82)  and  (86)  can  be  written  as  two  linear  algebraic 
systems  for  the  unknowns  (ptc)lCt jch  and 

m{K)px  +  At  ^  F(K,P)pk  —  F(P,K)pp  =  m(k')pkK,  for  all  K  €  fCi, .  (93) 

PzNiK) 

m(K)SK  +  &L  Y,  F(K'  P)6x  ~  K)$P  -  m(K)Hkn ,  for  all  K  €  K,h .  (94) 

P€Jf(K) 

If  a  volume  control  L  corresponds  to  a  node  which  is  on  the  boundary  then 
wo  set  pp  =  pQ,  Sp  —  0  and  we  eliminate  the  corresponding  equation. 
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Figure  9:  The  finite  element  mesh. 

3.5  Numerical  results 

In  this  section,  we  apply  the  material  and  computational  model  developed  to  simu¬ 
late  the  tunneling  phase  of  penetration  into  a  concrete  target  of  a  rigid  projectile. 
The  simulations  correspond  to  an  impact  velocity  of  V0  —  500  m/s  and  a  projectile 
of  mass  M  —  6  kg. 

3.5.1  Geometrical  and  physical  settings 

Preliminary  numerical  simulations,  in  which  a  pari  of  the  tunnel  was  considered  in 
V{)  have  shown  that  the  material  behind  the  projectile  is  rigid  and  thus  does  not 
affect  the  computations.  Hence,  to  reduce  the  computational  efforts,  the  control 
volume  Vo  is  chosen  such  that  it  does  not  touch  the  tunnel,  i.e.  F|  =  f). 

The  finite  element  mesh  of  the  domain  Z>0  in  cylindrical  coordinates  (in  the  Orz 
plane),  is  plotted  in  Figure  9.  It  has  33169  nodes  and  167666  triangular  elements. 
The  finite  volume  mesh  may  coincide  with  the  finite  element  triangulation. 

In  a  computation  of  the  entire  penetration  process,  the  initial  conditions  for  the 
tunnelling  phase  should  correspond  to  the  final  state  attained  in  the  previous  phase, 
called  "sabot  impact  phase”  (i.c.  the  phase  when  the  projectile  is  fully  embedded 
in  the  target  and  the  pusher  plate  hits  the  target  and  strips  off).  However,  in 
order  to  compute  this  state,  it  is  necessary  to  model  all  the  phases  of  penetration 
proccss( including  the  nose  penetration  phase  and  partially  restrained  penetration 
phase  corresponding  to  chipping  and  cratering  at  the  point  of  impact),  which  is 
beyond  the  scope  of  the  present  paper.  Thus,  in  this  work  we  consider  that  the 
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initial  conditions  v^p^ are  given  by  the  stationary  solution  corresponding  to  the 
striking  velocity  Vo  —  500  in/s.  The  solution  of  the  stationary  problem  is  computed 
using  the  techniques  developed  in  our  previous  work  [17], 

The  time  interval  for  the  computations,  [0,T],  is  chosen  such  that  at  /  —  7\  the 
velocity  of  the  projectile  is  :V(T)  =  Vmin  =  50  m/s.  This  choice  of  Vmm  is  related 
to  the  fact  that  for  very  low  velocities  (i,c.  V(l)  <  Vmm  =  50  m/s)  elastic  effects, 
which  were  neglected  in  our  model,  become  important  with  respect  to  the  plastic 
ones.  Note  however  that  the  phase  of  the  process  which  is  not  modeled  in  this  work 
(i.e.  for  V(t)  <  Vmin  =  50  m/s)  corresponds  to  less  than  5  %  of  the  entire  trajectory. 

The  specific  expressions  of  the  constitutive  functions  n(p,6)  and  pr[p)i  involved 
in  the  material  model  (equations  (G8)-(69))  are  taken  from  Holmquist  cl  ah  [30]. 
The  function  pr  is  expressed  in  terms  of  p  =  p/p°  —  1  as: 


* 


Pc(p)  -  P 


Pi'rujih 

perush 


Pc(p)  “  Pcruch  + 


[fl  flrrufth)  iptock  Pcruch) 


Ptock  Petits  h 


Pc(p)  =  A'!  j  -  +  A'a 

1  +  Plork 


P  ~  Plork 


f  "t"  Plork  . 


+  A;j 


P  ~  Plork 
“I-  Plork  . 


If  P  — 

if  P  €  ?  P<\  » 

if  n  >  ft,. 


(95) 

The  values  of  the  material  constants  involved  in  the  above  equation  fire:  p„ul,h  ~ 
1G  MPa,  ptock  =  800  MPa,  Pcru»h=  0.001,  plock  =  0.1,  //.  =  0.111.  Kj=  85  GPa, 
A'2=-171  GPa,  and  A';,=208  GPa  (see  [30]). 

The  yield  limit  in  shear,  n(p,8),  which  depends  on  both  the  compaction  level 
find  damage  is  expressed  as 


k(M)  =  ^o\/273nim{Sl„M,j4(l  -  5)  +  B{pc{p) / aa)s  } .  (96) 


with  cr0=  48  MPa,  Smai  =  7,  A  =  0.79,  B  =1.60,  and  N  =  0.61.  Following  [30],  it 
is  assumed  that  the  loss  of  cohesive  strength  of  the  concrete  material  is  related  to 
damage  from  irreversible  deformation,  which  evolves  according  to  the  following  law: 


as 

-  +  V.v*  = 


£>»[ 

Mp)  ' 


(97) 


where  A(/<)  —  y^2/3  inax{r0,  Dj[(pc(^)  +  T)fau\D* },  where  rih  Dt,  D2  are  constants 
(r0  =  0,01,  Z?)  =  0.04,  D-2  =  1  for  the  given  material)  and  T  is  the  maximum  tensile 
hydrostatic  pressure  that  the  material  can  withstand.  For  the  concrete  material 
studied.  T  =  4  MPa  ,  Note  that  using  the  continuity  equation  and  the  definition  of 
/i.  we  have: 

^  +  v  ■  Vji  =  -(1  +  p)  div(u), 

Thus,  the  model  for  concrete  is  of  the  general  form  (66)  with  n  =  2,  h  =  {p.S)  and 
H(D(v).h)  =  (-(1  +/i)div(o)  ,  |Z>'(T))|/A(/i)). 
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Figure  10:  Penetration  depth  as  a  function  of  time  for  different  time  steps. 

The  volume  and  shear  viscosity  coefficients  arc  considered  constant  and  have  the 
reference  values  A  =  0,002  MPa*s  and  tj—  0.001  MPa  s. 

FVictional  contact  is  assumed  on  the  boundary  P2  of  the  projectile  (sec  Figure 
1),  To  account  for  rate  effects  on  the  frictional  contact  between  the  concrete  target 
and  the  metallic  penetrator,  a  decaying  exponential  dependence  on  the  slip  rate 
proposed  in  [31]. 

■F(v)  =  Foo  +  exp(-cfv/vo){Fo  ~  *>)■ 

was  considered.  In  the  above  equation,  JF0  represents  the  static  value  of  the  friction 
coefficient,  while  is  its  dynamic  value,  while  Cj  is  a  constant.  Since  experimental 
determination  of  the  values  of  these  coefficients  is  very  difficult,  the  values  of  these 
coefficients  were  set  to:  J =  0  02,  —  0.2,  Vq  —  1  m/s  and  Cj  —  0.QL 

The  boundary  condition  un  —  0  prescribed  on  the  projectile  is  accurate  every¬ 
where  with  the  exception  of  two  small  zones.  The  first  one  is  located  on  the  nose  of 
the  projectile.  In  this  zone,  fracture  in  the  target  material  could  create  a  free  bound¬ 
ary,  A  second  zone,  which  is  located  on  the  shank  of  the  projectile,  corresponds  to 
the  separation  between  the  target  material  and  the  projectile.  The  formulation  and 
solution  of  these  two  free  boundary  problems  is  beyond  the  scope  of  this  study. 
However,  the  influence  of  these  two  zones  on  the  resistance  to  penetration  seems  to 
be  very  small. 

3,5.2  Time  step  sensitivity 

Let  first  check  the  time  step  sensitivity  of  the  proposed  explicit /implicit  discrctiza- 
tion  scheme.  To  this  end,  simulations  were  performed  using  time  steps  A/  varying 
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Figure1  11:  The  distribution  of  the  deformation  rate  |D(v)|  (  s~])  in  the  target* 
Left;  at  the  initial  stage  of  the  computations  (t  =  0,  V  —  500 m/s).  Right:  at  the 
final  stage  of  the  computations  (t  =  1,85ms,  V  =  50 m/s). 

from  5  ft  s  to  100  fi  s.  Figure  10  shows  the  evolution  of  the  penetration  depth  Dp  : 


versus  time  for  different  time  steps  A L  It  is  clear  that  even  for  a  coarse  time  step 
A /  —  50 fisy  which  corresponds  to  only  36  time  iterations,  the  penetration  depth  is 
accurately  predicted*  However,  for  time  steps  At  larger  than  100/is  the  time  dis¬ 
cretization  error  can  be  important*  We  conclude  that  the  time  discretization  scheme 
is  robust  and,  since  we  can  use  a  small  number  of  time  iterations,  the  computational 
cost  is  low. 

3.5.3  Fields  distribution  in  the  target 

In  Figure  11  is  shown  the  distribution  of  the  deformation  rate  |I>(t?)|  in  the  target 
corresponding  to  the  initial  and  final  stages  of  the  computations,  respectively  for  an 
impact  velocity  of  500  m/s.  It  is  seen  that  a  viscoplastic  zone  develops  around  the 
penetrator,  Outside  this  zone  there  is  no  deformation,  hence  the  target  is  rigid,  i.e. 
IDft?)!  vanishes.  It  can  be  observed  that  the  maximum  deformation  rate  is  reached 
on  the  nose  tip.  The  viscoplastic  zone  is  diminishing  during  the  penetration  process* 
Figure  12  (left)  shows  the  distribution  of  density  in  the  target  at  t  =  0.925ms,  tin 
intermediate  time  in  the  computation.  Note  the  existence  of  three  distinct  zones:  ( 1 ) 
mi  (almost)  incompressible  zone  around  the  penetrator  where  the  concrete  is  almost 
fully  compacted  {ft  >  /n)  ;  (2}  a  compressible  zone  where  the  density  is  greater  than 
p0?  the  density  of  the  intact  material*  mid  (3)  a  zone  that  remains  unaffected  by  the 
impact  event  where  the  density  is  everywhere  equal  to  It  is  worth  noting  that 
all  the  particles  ahead  of  the  projectile  which  are  at  a  distance  from  the  centerline 
less  than  the  projectile  radius  will  enter  the  fully  compacted  zone  (1), 

In  Figure  12  (right)  is  depicted  the  distribution  of  damage  (parameter  <$}  in  the 
target  at  mi  intermediate  time  of  the  computation:  t  —  0.925m Fully  damaged 
state  corresponds  to  &  =  1  while  S  =  0  corresponds  to  no-damage*  Let  note  that 
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Figure  12;  Left:  the  distribution  of  the  density  p  (  kg/m3)  in  the  target  at  t  — 
0.925 ms.  Right:  the  distribution  of  damage  6  in  the  target  at  t  =  0.925ms. 
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Figure  13:  The  distribution  of  the  rate  of  volumetric  deformation  div  v  (  s”1)  in  the 
target.  Left:  at  the  initial  stage  of  the  computations  (£  =  0,  V  =  500 m/5).  Right: 
at  the  final  stage  of  the  computations  (t  =  1,85ms,  V  =  50 m/s). 


according  to  the  model,  a  large  zone  of  the  target  is  fully  damaged.  The  extent 
of  the  fully  damaged  zone  varies  little  during  the  penetration  process  (i.e.  there  is 
little  variation  with  !/(£)). 

The  distribution  of  the  rate  of  volumetric  deformation  divt?  at  the  initial  and 
final  stages  of  the  computations  are  plotted  in  Figure  13.  We  note  that  around  the 
penetrate r  an  almost  fully  compacted  state  is  achieved,  the  material  being  almost 
incompressible.  The  maximum  compaction  is  reached  in  the  nose  zone,  ahead  of  the 
projectile. 

Let  us  analyze  now  the  stress  distribution  in  the  target  and  determine  the  zones 
where  tensile  failure  may  occur.  It  is  assumed  that  concrete  tensile  failure  can  be 
modeled  with  the  classic  maximum  tensile  strength  criterion  (see  [32]),  i.e,  the  ma¬ 
terial  fails  if  the  maximum  of  the  principal  stresses  reaches  a  critical  limit  denoted 
of.  From  the  data  reported  in  [30]  oy  =  4  MPa.  From  the  analysis  of  the  stress 
distribution  in  the  target,  we  found  that  the  only  principal  stress  which  can  have 
positive  values  (i.e.  is  tensile)  is  oe&  (sec  regions  in  green  in  figure  14).  As  shown 
in  figure  14  tensile  failure  occurs  in  a  wide  region  (colored  in  red)  ahead  of  the 
projectile.  This  region  is  smaller  for  low  projectile  velocities  i,e,  at  the  end  of  the 
penetration  process.  From  this  analysis,  fundamental  understanding  of  the  stability 
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Figure  14:  The  distribution  of  ae&  in  the  target.  In  green  the  region  which  are  in 
tensile  (cT00  >  0)  and  in  red  the  regions  where  tensile  failure  occurs  (am  >  aj). 
Left:  at  the  initial  stage  (t  —  0,  V  —  500m fs).  Right:  at  the  final  stage  of  the 
computations  (t  —  1.85ms,  V  =  50 m/s). 

of  the  pcnct rater  trajectory  can  be  gained.  Indeed,  since  the  tensile  failure  could 
occur  only  along  Ozr  planes,  orthogonal  to  which  are  not  symmetric  with  re¬ 
spect  to  the  pcnctrator  centerline,  anisotropy  is  generated  in  the  target.  Trajectory 
deviations,  observed  at  high  penetration  velocities,  may  be  related  to  this  induced 
anisotropy. 

3*5,4  Parameter  sensitivity  of  the  penetration  depth 

As  already  mentioned,  one  of  the  many  challenges  associated  with  modeling  penetra¬ 
tion  is  that  only  very  limited  experimental  information  such  as  projectile  trajectory 
and  deceleration  data  can  be  gathered  during  penetration.  To  date,  there  is  limited 
data  from  which  the  yield  strengths  may  be  inferred  and  thus  the  cut-off  yield  limit 
Smax  involved  in  the  expression  (96)  of  the  shear  yield  limit  k  cannot  be  determined 
explicitly.  Likewise,  it  is  very  difficult  to  evaluate  experimentally  the  viscous  coef¬ 
ficients  A  and  rp  and  the  coefficient  cj ,  involved  in  the  friction  law.  Because  of  the 
lack  of  test  data,  a  parametric  study  was  conducted  to  estimate  the  sensitivity  of 
the  depth  of  penetration  predictions  to  these  constants. 

First  it  was  assessed  the  importance  of  plastic  versus  viscous  effects  oil  the  pen¬ 
etration  depth*  To  better  analyze  these  effects,  friction  was  completely  neglected, 
i.c*  we  set  T  =  0. 

To  estimate  the  influence  of  viscosity  on  the  penetration  depth,  we  have  set  all 
the  parameters  involved  in  the  expression  of  k ,  the  yield  limit  in  shear,  anti  of  pr 
to  be  equal  to  the  reference  values  and  we  have  performed  a  scries  of  numerical 
simulations  in  which  we  varied  the  values  of  the  viscous  parameters  t)  and  A.  In 
Figure  15  are  plotted  the  evolution  of  the  penetration  depth  t  — *  Z?p(f)  with  time 
for  different  values  of  the  shear  viscosities  tj.  In  all  these  calculations,  the  volume 
viscosity  A  was  held  constant  (i.c*  equal  to  the  reference  value).  It  is  clearly  seen 
that  for  r;  less  than  1000  Pa  s,  the  plastic  properties  of  the  target  are  dominant  while 
the  viscous  effects  are  negligible. 


50 


tin*  (m« » 


Figure  15:  Evolution  of  penetration  depth  t  — *  Dp(t)  for  different  shear  viscosities 
*]• 


However,  for  r?  greater  than  1000  Pa^s,  the  viscous  effects  are  more  pronounced. 
Indeed,  in  the  range  of  103  —  10s  Pa-s  a  variation  of  two  orders  of  magnitude  of  the 
value  of  rj  gives  rise  to  an  order  of  magnitude  variation  in  penetration  depth  . 

We  have  also  computed  the  evolution  of  the  penetration  depth  t  —*  Dp(t.)  for 
different,  values  of  the  volume  viscosity  A  in  the  range  of  102  —  105  Pa-s.  In  all  the 
calculations,  the  shear  viscosity  tj  was  held  constant  (equal  to  a  reference  value). 
We  have  found  that  the  volume  viscosity  has  no  influence  on  the  penetration  depth. 
An  explanation  of  this  insensitivity  to  the  volume  viscosity  could  be  the  fact  that 
around  the  projectile  the  target  material  is  almost  incompressible  (see  figures  12 
and  13),  hence  the  volume  strain  rate  div(u)  is  very  small.  That  means  that  the 
stress  distribution  around  the  projectile  is  no  affected  by  the  volume  viscosity  A 
which  multiplies  div(r).  Since  the  resistance  force,  involved  in  the  calculation  of  the 
penetration  depth,  is  computed  based  only  on  the  distribution  of  the  stress  tensor 
<t  on  the  projectile  surface,  we  can  deduce  that  one  cannot  expect  any  variation  of 
the  resistance  force  with  respect  to  the  volumetric  viscosity  A. 

In  addition,  since  it  is  very  difficult  to  determine  experimentally  the  viscosity 
coefficients  of  concrete  for  dynamic  loading  (high  rate  of  deformations  and  high 
stresses),  the  sensitivity  analysis  performed  is  important  in  selecting  the  numerical 
values  of  these  constants  that  will  provide  the  best  fit  to  the  penetration  depths 
which  arc  directly  measured  in  penetration  tests. 

In  order  to  sec  the  effects  of  the  plastic  parameters,  we  set  the  the  viscous  param¬ 
eters  A  and  t)  equal  to  the  reference  values  and  wc  performed  a  scries  of  numerical 
simulations  for  different  values  of  the  cut-off  yield  limit  Smax,  parameter  involved 
in  the  expression  of  the  yield  limit  k  which  cannot  be  obtained  from  experiments  in 
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Figure  16:  Evolution  of  penetration  depth  t  — ►  Dp(t)  for  different  values  of  the 
cut-off  yield  limit  Smax. 

tin  explicit  manner,  As  in  the  previous  calculations,  we  have  completely  neglected 
friction  he.  we  took  {F  =  0). 

In  Figure  16  is  plotted  the  evolution  of  the  penetration  depth  /  Dp(t)  for 
different  values  of  the  cut-off  yield  limit  SmtLt.  Note  that  for  greater  than 
9,  there  is  little  influence  on  the  penetration  depth.  That  means  that  near  to  the 
projectile,  the  pressures  D{pc{p)/a^)N  tire  less  than  the  cut-off  value  S1UUX.  If  5„iax 
is  less  than  9,  then  Smta  constitutes  a  limit  on  the  strength  that  can  be  developed 
in  the  material  (i,c,  really  works  as  a  cut-off  value  for  the  yield  strength).  In  this 
case,  it  was  found  that  the  penetration  depth  is  very  sensitive  to  the  yield  stress. 
Indeed,  if  Smax  is  three  times  smaller  then  the  final  penetration  depth  is  three  times 
larger. 

Next,  the  effect  on  the  penetration  depth  of  the  coefficient  C/T  involved  in  the 
slip  rate  weakening  friction  law  was  analyzed.  To  this  end,  we  held  all  the  other 
parameters  constant  (he.  we  set  JF0  and  and  all  the  coefficients  involved  in  the 
model  for  the  target  equal  to  reference  values).  In  Figure  17  is  plotted  the  evolution 
of  penetration  depth  t  —>  Dp(t)  for  different  values  of  C/.  We  note  that  if  c.j  is 
greater  than  0.Q1  it  doesn't  have  any  effect  on  the  final  penetration  depth.  In  this 
case,  the  friction  coefficient  reaches  its  dynamical  value  almost  at  every  point 
on  the  projectile.  If  0/  is  less  than  0.01,  then  the  slip  rate  weakening  of  the  friction 
coefficient  becomes  apparent,  i.c.  the  penetration  depth  is  also  sensitive  to  the  value 
of  the  static  friction  coefficient 
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Figure  17:  Evolution  of  penetration  depth  t  — *  Dp(t)  for  different  values  of  the 
exponential  decrease  coefficient  cj,  involved  in  the  slip  rate  weakening  friction  law. 

4  Conclusions 

We  have  developed  computational  methods  for  modeling  the  penetration  of  a  rigid 
projectile  into  porous  geological  and  cementitious  targets.  To  capture  the  solid-fluid 
transition  in  behavior  at  high  strain  rates  observed  in  such  media  and  account  for 
damage/plasticity  couplings  and  how  these  dissipative  mechanisms  arc  influenced 
by  the  strain  rate,  a  compressible  rigid  viscoplastic  fluid  constitutive  equation  was 
considered.  A  hybrid  time-discretization  was  used  to  model  the  non-station  ary  flow 
of  the  target,  i.e.  an  explicit  Euler  method  for  the  projectile  equation  and  a  forward 
(implicit)  method  for  the  target  boundary  value  problem.  At  each  time  step,  a  mixed 
finite-element  and  finite-volume  strategy  was  used  to  solve  the  "target”  boundary 
value  problem. 

Specifically,  variational  formulations  and  algorithms  for  solving  the  minimiza¬ 
tion  problems  with  non-linear  constraints  for  the  velocity  field  via  the  finite  element 
approximation  were  developed  while  finite-volume  techniques  were  adopted  for  solv¬ 
ing  the  hyperbolic  mass  conservation  and  damage  evolution  equations.  Since  in 
the  constitutive  description  of  the  target  material  a  rigid  unloading  hypothesis  was 
considered,  the  variational  inequality  to  be  solved  has  unilateral  constraints  for  the 
velocity  divergence.  Additional  difficulties  were  related  to  ensuring  incompressibility 
in  the  regions  where  the  material  is  unloading.  Also,  due  to  plastic  and  frictional 
contributions,  the  resulting  functional  is  not  smooth.  To  overcome  these  difficulties, 
the  decomposition-coordination  formulation  coupled  with  the  augmented  lagrangian 
method  proposed  for  the  incompressible  Bingham  fluid  by  [19,  20]  was  adapted  for 
the  material  model  considered  in  this  work. 
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Numerical  simulations  of  penetration  into  concrete  were  performed.  It  was  shown 
that  using  the  model  and  numerical  techniques  presented  in  this  paper,  it  is  possible 
to  predict  the  density  changes  around  the  penetration  tunnel,  the  shape  and  location 
of  the  rigid/plastic  boundary,  the  extent  of  damage  due  to  air  void  collapse  in 
such  materials.  According  to  our  calculations,  tensile  failure  occurs  ahead  of  the 
penctrator  and  along  planes  which  arc  not  symmetric  with  respect  to  the  pcnctrator 
centerline.  This  induced  anisotropy  may  explain  trajectory  deviations. 

One  of  the  many  challenges  associated  with  simulating  penetration  in  cemen¬ 
titious  materials  is  related  to  the  lack  of  experimental  data  for  the  high  pressures 
and  high  strain  rates  involved.  Thus,  some  of  the  constants  involved  in  the  model 
cannot  be  determined  directly  from  experimental  data.  A  parametric  study  was  per¬ 
formed  to  asses  the  sensitivity  of  the  predicted  penetration  depth  to  these1  constants 
(i.e.  to  shear  and  volumetric  viscosities,  cut-off  yield  limit,  and  dynamic  friction 
coefficients).  The  results  of  this  parametric  study  provides  insights  into  the  relative 
importance  of  the  plastic  and  viscous  effects  on  the  penetration  process.  Finally,  by 
conducting  a  time  step  sensitivity  study  it  was  shown  that  the  numerical  model  is 
robust  and  computationally  inexpensive. 
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